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CHAPTER  1 

TNTRQDUCTION 

Human  beings  have  been  trying  to  do  things  in  the  most  efficient 
way  possible  since  the  beginning  of  civilization.  From  the  pre- 
historic period  to  the  present  age  the  accomplishment  of  tasks  with 
the  least  expenditure  of  energy  or  time,  has  been  a  motivating  factor 
for  expanding  the  technological  base. 

At  the  beginning  of  the  industrialized  age,  human  labor  became  an 
integral  part  of  what  was  then  "the  art  of  manufacturing."  A  few 
people  were  able  to  produce  products  at  a  higher  rate  than  others  and 
they  were  called  "experts."  F.W.Taylor  [1],  in  the  early  1900' s, 
realized  that  there  were  certain  motions  used  by  these  experts  that 
made  them  perform  their  tasks  faster  than  others.  Taylor  fashioned 
and  promoted  the  idea  of  making  others  follow  similar  motions  to 
optimize  the  production  rate  and,  hence,  developed  the  science  of  work 
study,  which  played  an  important  role  in  the  optimization  of 
production  in  flow  line  systems  of  manufacturing. 

As  technology  changed  and  flow  lines  gave  way  to  fully  automated 
lines,  the  need  to  use  the  most  efficient  movements  of  machines  has 
become  a  primary  concern.  The  branch  of  mathematics  known  as 
"optimization"   has   found  many  new  applications  in  tuning  the 


production  processes.  Optimization  mathematics  facilitated  design  of 
manufacturing  machinery  which  operated  at  increasingly  higher  rates 
and  efficiencies.  However,  a  drawback  of  this  fixed  manufacturing  was 
that  in  order  to  manufacture  a  different  product,  the  entire  assembly 
line  had  to  be  re -tooled  to  accomodate  the  new  product.  With  the 
introduction  of  computers  to  the  manufacturing  environment  it  became 
possible  to  route  different  products  through  the  same  line  with 
minimal  changes.  Although  computer  control  made  the  transition  from 
one  production  schedule  to  another  relatively  easy,  the  need  to 
optimize  the  production  of  a  single  product  still  remained. 

There  is  a  cost  of  manufacturing  associated  with  any  production 
schedule.   The  first  task  in  optimizing  or  minimizing  cost  is  to 
express  cost  as  a  function  of  those  quantities  which  we  are  free  to 
vary  and  can  control.   The  two  most  accessable  of  these  quantities  are 
energy  and  time  consumption.  We  can  minimize  consumption  of  energy  by 
reducing  the  rate  of  production  or  by  selecting  more  efficient  drives. 
Obviously,   reducing  production  rates  too  greatly  has  the  undesirable 
consequence  of  curtailing  output.   Minimizing  the  consumption  of  time 
can  be  achieved  by  minimizing  the  cycle  time  of  the  individual 
components  of  the  manufacturing  system  as  discussed  above.   Time  of 
production,  subject  to  the  available  resources,  provides  a  viable  area 
in  which  to  achieve  lower  costs.   This  can  be  done  by  optimizing  the 
cycle  time  of  the  individual  components  of  the  manufacturing  system. 

One  avenue  through  which  the  optimization  or  in  this  case 
reduction  of  cycle   time  can  be  achieved  is  to  increase  the  speed  at 


which  the  production  machinery  or,  in  the  context  of  this  work, 
robots  operate.  The  problem  of  cycle  time  optimization  of  mechanical 
manipulators  falls  into  the  general  category  of  minimum  time  control 
problems  with  or  without  path  constraints.  Parenthetically, 
production  operations  have  path  constraints  in  order  to  avoid 
collisions  with  the  assembled  part  or  other  obstacles  in  the  work 
space.  On  the  other  hand  loading  and  unloading  tasks  in  metal  working 
operations  and  missile  interception  problems  are  examples  of  tasks 
without  path  constraints.  Extensive  work  has  been  done  in  determining 
the  path  constrained,  minimum  time  motion  [2]  -  [5].  This  work  will 
deal  with  the  non-path- constrained  minimum  time  problem. 

Time  optimal  control  without  path  constraints  can  be  defined  as 
the  control  function  which  forces  a  system,  described  by  a  system 
equation  and  having  bounded  control,  from  one  point  to  another  in 
minimum  time.  The  system  equation  is  normally  a  differential 
equation.  Systems  which  are  defined  by  linear  differential  equations 
are  called  linear  systems  and  those  defined  by  nonlinear  diffential 
equations  are  called  nonlinear  systems.  The  early  development  of  time 
optimal  control  problems  without  path  constraints  was  done  almost 
exclusively  in  terms  of  linear  systems.  Nonlinear  systems  were  later 
solved  by  using  some  of  the  techniques  developed  in  the  solution  for 
linear  systems.  Some  of  the  solution  techniques  developed  for  linear 
and  nonlinear  systems  will  be  presented  here.  It  should  be  noted  that 
even  though  the  theory  might  look  simple,  in  practice  it  is  difficult 
to  apply  these  solution  techniques  to  all  systems. 


LITERATURE  SURVEY 

T.TNEAR  SYSTEMS 

The  solution  techniques  '  developed  for  linear  systems  can  be 

broadly  classified  as  analytical  and  numerical. 

ANALYTICAL  MODELS 

The  early  development  of  the  solution  techniques  for  linear 
systems  involved  finding  the  analytical  solution.  Significant  among 
the  early  analytical  models  was  that  presented  by  Athanassiades  and 
Smith  [6].  Athanassiades  and  Smith  developed  a  minimum  time  control 
system  for  an  Nth  order  linear  system  with  negative  real  poles.  They 
transformed  the  state  equations  into  a  diagonalized  form,  solved  for 
time  as  a  logarithmic  function  and  obtained  the  switching 
hypersurface.  Based  on  the  position  of  the  state  point  relative  to 
the  switching  surface,  they  were  able  to  determine  the  control.  The 
control  they  determined  would  drive  the  system  from  the  initial  state 
to  the  phase  space  origin  in  the  least  time  possible  given  the 
capacity  of  the  available  drives.  At  about  the  same  time,  working 
independently,  Desoer  and  Wing  [7]  came  up  with  a  similar  technique 
for  discrete  systems. 

Later,  Ryan  [8]  developed  algebraic  expressions  for  the  minimum- 
time  isochronal  surfaces  for  third  order  systems  with  real  eigen 
values  and  a  single  saturating  input.  He  integrated  the  system 
equations  using  simplified  controls  and  solved  for  the  final  tii 
This  becomes  difficult  for  higher  order  and  nonlinear  systems. 
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Another  interesting  approach  was  presented  by  Fuller  [9].  He 
computed  the  time-optimal  control  required  for  points  where  all  but 
one  coordinate  are  zero  by  backward  tracing  of  the  trajectories  for  a 
linear  plant  with  distinct  real  eigenvalues  and  a  scalar  saturating 
control.  The  trajectory  was  divided  into  m  time  intervals  at  which 
the  control  changes  sign.  Backward  integration  of  the  system 
equations  from  the  origin  then  provides  the  optimal  control.  He 
suggested  that  for  state  points  away  from  the  axes  the  control  can  be 
determined  from  the  dominant  state  coordinate.  However  there  is  a 
possibility  of  chatter  near  the  switching  curve  in  this  case. 

Sebakhy  and  Abdel-Moneim  [10]  and  Rao  and  Janakiram  [11] 
presented  algorithms  to  compute  the  time-optimal  control  laws  to  drive 
closed  loop  discrete  linear  systems  to  the  origin  of  the  phase  space. 
Sebaky  and  Abdel-Moneim  defined  the  system  by 

x(t+l)  -  A  x(t)  +b  u(t)  (1-1) 

u(t)  -  F  x(t)  (1-2) 

where   x  is  an  n  dimensional  state  vector,  u  is  the  control  vector, 

and  A,  b  and  F  are  real  matrices.   The  system  was  assumed  to  be 

controllable.   They  computed  the  feedback  law  F  from 

(A  +  b  F)'^  -  0  (1-3) 

such  that  r  is  minimum. 

Rao  and  Janakiraman  presented  a  solution  technique  for  discrete 
systems  specified  in  the  z  -  domain  of  the  form, 


Uz)_.^J^ ;p<n  (1.4) 

M(z)    n      . 

I     \   ^' 
i-0  ^ 

where  X(z)  and  M(z)  are  the  z  -  transforms  of  the  output  and  input 
respectively.  They  defined  the  state  variables  as  the  output  at  n 
consecutive  time  intervals  and  solved  for  the  input  as 

[&]    [ffl]  -  -[b]  [x(0)]  (1.5) 

giving, 

[m]  -  -[a"^]  [b]  [x(0)  ]  (1.6) 

where  [a]  and  [b]  were  obtained  by  converting  the  system  to  its 
difference  equations.  In  cases  where  there  was  a  saturation  limit  on 
the  controls,  the  input  was  set  at  that  limit. 

NUMERICAL  MODELS 

With  the  advent  of  the  computer,  numerical  solution  techniques 
became  popular.  The  solution  methodologies  for  time -optimal  control 
problems  were  directed  towards  iterative  techniques. 

Knudsen  [12]  developed  an  iterative  procedure  for  computing  the 
time -optimal  controls  for  the  generalized  state  equations  of  the  form 

i  -  A  X  +  b  u  (1.7) 

where  A  and  b  are  constant  matrices,  x  the  state  vector,  x  the  rate  of 

change  of  the  state  vector,  and  u  is  the  control  driving  the  system. 
The  control  satisfies  the  inequality  constraint  of 

|u|  <  1.  (1-8) 


He  parameterized  the  minimum  time  control  in  terms  of  the  initial 
conditions  of  the  system's  adjoint  equations.  Knudsen  derived  a 
function  relating  the  initial  state  of  the  system  to  the  adjoint 
system's  initial  condition  and  to  the  time  required  by  the  optimal 
control  to  drive  the  state  to  zero.  This  function  was  used  to  solve 
the  problem  iteratively.  Nagata,  Kodama  and  Kumagai  [13]  developed  an 
iterative  procedure  to  solve  discrete  state  variable  formulations. 

Lastman  [14]  developed  a  shooting  method  for  two  point  boundary 
value  problems  arising  from  bang-bang  control  problems .  He  guessed 
the  final  time  and  the  initial  values  of  the  multipliers  and 
iteratively  improved  the  solution  by  integrating  forward  in  time.  He 
also  showed  that  nonlinear  systems  can  be  solved  by  this  method. 

Another  interesting  numerical  technique  developed  by  Larson  [15] 
and  Yastreboff  [16]  consisted  of  iteratively  changing  the  time 
interval  between  switchings.  Larson  assumed  arbitrary  switching 
intervals  and  controls  and  determined  an  initial  trajectory.  He  then 
determined  a  sequence  of  trajectories  by  correcting  the  previous 
trajectory  so  that  the  respective  sequence  of  end  points  of  the 
trajectories  will  approach  the  desired  point  in  state  space.  The 
correction  routine  used  a  first  order  Taylor  series  approximation 
about  the  present  switching  interval  times  and  assured  that  each 
trajectory  was  time  optimal  from  the  initial  to  whatever  final  point 
it  reached.  The  corrections  were  proportional  to  the  distance  by 
which  the  final  state  was  missed.  It  should  be  noted  that  the  terms 
neglected  in  the  Taylor  series  expansion  are  not  necessarily  small. 


The  optimal  control  was  obtained  when  the  correction  routine  brought 
the  end  point  of  the  trajectory  near  the  desired  end  point  in  state 
space.  The  method  was  applied  to  linear  systems  with  complex  eigen 
values  and  to  time  varying  systems. 

Yastreboff  [16]  formulated  his  time  interval  adjustment 
technique  for  linear  plants  with  constrained  control  amplitudes.  He 
arbitrarily  chose  n  switching  times  and  a  control  function  which 
brought  the  plant  to  the  desired  terminal  state.  The  switching  times 
and  control  are  then  adjusted  so  that  the  controls  approximate  a  bang- 
bang  form.  He  adjusted  the  control  by  a  linearization  technique  and  a 
logarithmic  technique.  The  linearization  technique  was  actually  a 
variational  technique.  He  took  the  first  variation  of  the  system 
equation  integrated  forward  in  time  and  set  the  variation  depending 
upon  the  final  state  to  zero.  The  problem  was  then  solved 
iteratively..  The  logarithhmic  technique  consisted  of  approximating 
the  variations  with  logarithmic  functions.  The  logarithmic  approach 
converged  faster,  but  determining  this  function  for  other  systems 
would  be  a  difficult  task. 

NONLINEAR  SYSTEMS 

Most  physical  systems  are  nonlinear  in  nature.  The  solution 
techniques  developed  for  linear  systems  are  useful  in  that  they  can  be 
used  to  solve  nonlinear  systems  by  linearized  approximations  in  the 
case  of  analytical  solutions  or  by  extending  them  directly  in  the  case 
of  numerical  techniques  as  in  Lastman  [14]. 


LINEARIZED  APPROXIMATIONS 

Notable  among  the  linearized  aproximation  techniques  used  was 
that  developed  by  Kahn  and  Roth  [17].   They  investigated  a  three 
degree  of  freedom  manipulator  and  obtained  a  suboptimal  solution  to 
the  minimum  time  problem  by  using  a  linearized  set  of  the  equations  of 
motion.   The  equations  were  decoupled  using  suitable  transformations. 
Approximations  were  made  to  the  gravity  and  velocity  terms  to  reduce 
the  problem  to  double  integrator  formulations  on  each  axis.   The 
method  gave  reasonably  good  results  on  those  axes  which  were  loosely 
coupled  with  others,  but  on  the  axes  which  had  a  higher  degree  of 
coupling,  the  error  between  the  actual  solution  and  the  solution  given 
by  the  linearized  equations  were  as  high  as  68%.   It  should,  however, 
be  noted  that  the  state  equations  were  simultaneous  coupled  nonlinear 
differential  equations. 

Wen  and  Desrochers  [18]  also  presented  two  linearized  solution 
techniques,  namely,  the  method  of  averaging  dynamics  and  the  method 
of  linear  equivalence.  The  method  of  averaging  dynamics  assumed  that 
the  non- linear  structure  was  constant  in  time  and  used  a  double 
integrator  solution  based  on  this  structure.  However,  velocity 
constraints  could  not  be  enforced.  The  method  of  linear  equivalence 
assumed  that  the  non-linear  structure  was  known  and  the  system  was 
decoupled  into  double  integrators.  The  torque  was  calculated  by  using 
the  non-linear  structure. 


NUMERICAL  TECHNIQUES 

Analytical  solutions  for  nonlinear  systems  are  very  difficult  to 
obtain.  Numerical  techniques  are  easier  to  apply,  however, 
computational  time  is  a  major  concern  in  this  area. 

Kahn  and  Roth  [17]  presented  a  numerical  solution  technique  to 
the  complete  nonlinear  problem  to  verify  their  linearized  equations. 
The  method  involved  making  a  guess  on  the  unknown  values  of  the 
adjoint  system  variables  at  the  final  time  and  integrating  both  the 
state  and  the  costate  equations  back  in  time  until  the  initial  state 
was  reached.  Iteratively  changing  this  guess  on  the  final  value  of 
the  costate  variables  produced  the  exact  control.  This,  however, 
cannot  be  used  for  an  on-line  control  scheme  which  is  required  in  most 
modern  production  systems  due  to  the  length  of  the  calculation. 

Shetty  [19]  presented  a  finite  element  approach  to  the  minimum 
time  problem  of  a  two  degree  of  freedom  robot.  He  formulated  a  finite 
element  model  of  the  optimization  equations  obtained  from  the 
Hamiltonian.  The  position  and  velocities  of  the  two  axes,  the 
multipliers  and  their  first  derivatives  and  the  elemental  time  were 
used  as  unknown  variables.  He  used  a  grid  search  method  to  get  a 
reasonably  good  guess  on  the  final  time  and  used  this  in  the  finite 
element  model.  The  nodal  variables  and  the  elemental  time  were 
iteratively  improved.  He  used  a  continuous  time  method  similar  to  the 
technique  presented  by  Lastman  [14]  to  verify  his  results. 

Davison  and  Monro   [20]   and  Wen  and  Desrochers  [21]  presented 
time  interval  optimization  techniques  for  nonlinear  systems.   Davison 
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and  Monro  presented  a  computational  technique  for  finding  the  time- 
optimal  controls  of  non-linear  time  varying  systems  with  single 
inputs.  They  assumed  the  number  of  switchings  and  the  time  of  each 
interval  and  formulated  a  composite  criterion  function,  which  was 
minimized  to  get  the  switching  times.  They  assumed  an  arbitrary 
number  of  switching  intervals  N  for  the  control  input  u(t)  and  made 
initial  guesses  on  the  switching  intervals  T^^'s  and  the  initial 
control.  The  response  was  calculated  based  on  these  values  and  the 
values  of  N  and  T's  were  refined  using  Rosenbrock's  Method  [22]  and  a 
composite  cost  function  of  the  form, 

J  -  c  T  +  x'^(t)  x(t)  (1-9) 

where  c  is  some  weighting  factor  (usually  c-1)  and 

N 
T-  y  T,  .  (1-10) 

i-1   ^ 

The  time  interval  optimization  technique  presented  by  Wen  and 
Desrochers  [21]  was  developed  for  nonlinear  systems  with  multiple 
inputs.  This  is  significant  due  to  the  fact  that  references  [15], 
[16],  and  [20]  considered  only  systems  with  single  inputs.  They 
guessed  the  N  switching  times  (t^,  ^^,  .  ..  t^)  where  a  switching  time 
is  defined  as  the  time  when  any  one  element  of  the  control  vector 
switches.  They  also  guessed  the  initial  control  vector,  the  final 
time  and  the  order  in  which  the  controls  switch.  They  iteratively 
improved  the  solutions  by  change  in  the  final  state  due  to  a  change  in 
an  earlier  switching  time.   The  method  worked  well  when  the  number  of 
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switches  assumed  is  greater  than  the  actual  number  and  with  some 
correct  intuitive  guesses. 

NOVEL  METHODS 

A  few  novel  methods  have  also  been  presented  to  reduce  the 
computational  effort  of  the  minimum- time  problem.  Goor  [23]  converted 
a  three  degree  of  freedom  robot  problem  into  three  separate  problems 
and  defined  the  path  based  on  the  slowest  axis.  The  problem  was 
converted  to  one  in  which  the  'jerk',  which  was  caused  by  the 
switching  of  the  controls,  was  bounded  to  reduce  wear  and  tear.  The 
solution  to  this  simplified  system  was  fovind  for  each  of  the  axis  and 
the  control  was  defined  based  on  the  slowest  axis.  This  does  not  give 
us  the  actual  minimum  time  control  or  use  the  non- linear  dynamics  to 
advantage . 

Luh  and  Shafran  [24]  presented  another  interesting  approach. 
They  used  two  succesive  least-squares-fit  approximations  to  get  the 
isochrones  of  the  system.  The  minimum-time  isochrones  were  computed 
in  a  discrete  region  by  using  some  known  method.  The  isochrones  were 
approximated  by  a  hyperellipsoidal  surface  in  terms  of  the  minimal 
time  and  the  state  variables.  The  coefficients  of  the 
hyperelliposoidal  surfaces  were  then  ordered  according  to  their 
corresponding  minimal  times  and  they  were  approximated  as  a  set  of 
continuous  functions  of  time.  These  approximate  functions  were  used 
to  generate  the  controls.  The  first  two  least-squares-fits  were  done 
off  line  and  only  the  approximate  function  was  used  for  real  time 
computation  purposes.   The  accuracy  of  the  isochronal  hypersurf aces , 
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however,   depends  on  the  number  of  state  points  chosen  for  the  initial 
leas t -  squares -  fit . 

PRESENT  APPROACH 
The  solution  techniques  available  for  time  optimal  control 
problems  without  path  constraints  have  been  discussed.  However,  there 
is  one  drawback  or  another  in  most  of  these  methods.  Analytical 
solutions  are  extremely  difficult  to  obtain  except  for  very  simple 
systems.  The  numerical  computational  techniques  work  well  with  linear 
systems  where  the  state  transition  matrix  is  known.  Even  then 
convergence  to  the  absolute  minimum  requires  intuitive  guessing.  In 
the  case  of  nonlinear  systems,  the  state  transition  matrix  is  not 
available,  the  integration  interval  and  stability  become  greater 
factors,  and  intuitive  guessing  becomes  extremely  difficult.  These 
factors  makes  numerical  techniques  more  difficult  to  apply.  Sometimes 
it  takes  an  enormous  amount  of  time  and  effort  to  come  up  with  the 
optimal  control  and  hence  it  can  only  be  used  in  path  planning. 
Another  important  factor  to  be  considered  is  that  the  computations 
have  to  be  repeated  for  each  initial  position. 

Due  to  the  problem  of  computational  time  and  repeated 
computations  required  for  various  initial  positions,  it  is  necessary 
to  come  up  with  an  efficient  method  for  an  online  control  scheme. 
Analysis  of  the  problem  in  the  state  space  domain  will  reduce 
unnecessary  computations  for  changes  in  initial  positions.  Luh  and 
Shafran  [24]  presented  one  such  approach.  However,  as  stated  before, 
the  accuracy  of  their  analysis  depended  upon  the  number  of  state 
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points  chosen.  The  purpose  of  this  work  is  to  provide  an  efficient 
method  to  solve  the  complete  minimum  time  problem  for  real  time 
applications  and  to  examine  the  feasibility  of  a  continuum  approach 
towards  such  a  solution. 

The  continuum  approach  presented  here  treats  the  state  variables 
as  independent  variables  in  the  state  space.  The  final  time  is  then 
evaluated  as  a  function  of  the  state  variables.  This  choice  of 
variables  enables  us  to  treat  both  linear  and  nonlinear  systems  in 
essentially  the  same  way.  The  continuum  relations  have  been  derived 
and  an  approximate  solution  has  been  obtained  for  a  double  integrator 
problem  using  a  finite  element  analysis. 


14 


CHAPTER  2 

THE  TIME  OPTIMAL  CONTROL  PROBLEM 
OPTIMAL  CONTROL  THEORY 

Classical  control  systems  design  is  based  on  acceptable 
performances  defined  in  terms  of  time  and  frequency  domain  criteria 
such  as  rise  time,  settling  time,  peak  overshoot,  gain,  and  phase 
margin,  and  bandwidth.  Modern  technology  demands  complex  multiple - 
input,  multiple -output  systems  with  radically  different  performance 
criteria.  With  the  development  of  the  digital  computer,  optimal 
control  theory  has  been  used  in  the  design  of  these  systems. 

The  objective  of  optimal  control  theory  [25]  is  "to  determine  the 
control  signals  that  will  cause  a  process  to  satisfy  the  physical 
constraints  and  at  the  same  time  minimize  (or  maximize)  some 
performance  criterion." 

In  order  to  evaluate  the  performance  criterion  (or  index)  and  to 
determine  the  optimal  control,  the  designer  must  have  a  complete 
knowledge  of  the  mathematical  description  (or  model)  of  the  process  to 
be  controlled,  a  statement  of  the  physical  constraints,  and  a 
specification  of  the  performance  index. 
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THE  MATHEMATICAL  MODEL 

The  mathematical  modelling  of  the  system  is  an  important  aspect 

of  any  control  problem.   The  model  should  be  an  accurate,  but  simple 

description  of  the  physical  system  and  should  predict  the  response  to 

all  anticipated  inputs  within  reasonable  accuracy.   If  the  states  of 

the  system  at  any  time,  t,  are  Xj^(t) ,  X2(t),  x^Ct) x^(t)  and 

the  r-nntrol  inputs  to  the  system  at  any  time,  t,  are  u^(t),  U2(t),  .  . 
,u  (t),   the  system  may  be  described  by  n  first  order  differential 
equations  of  the  form 

x^(t)  -  a^(x^(t),  X2(t) '^n^^^' 

Uj^(t),  U2(t),  ....  Uj^(t).  t), 

X2(t)  -  a2(x^(t),  X2(t) ^n^'^^ ' 

u^(t),  U2(t) Uj^(t),  t), 

(2.1) 


3c^(t)  -  a^(x^(t),  X2(t) '^n^^^' 

Uj^(t),  U2(t),  ....  u^(t),  t) 


Then,  let 
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x(t) 


Xj^Ct) 


X2(t) 


x^(t) 


be  defined  as  the  state  vector,  and  let 


Uj^Ct) 


(2.2) 


u(t) 


u, 


2(t) 


(2.3) 


be  defined  as  the  rnnttrol  vector.   The  system  can  now  be  defined  in 
the  vector  form  as 


x(t)  -  f(x(t),  u(t),  t), 


(2.4) 


where  f  is  a  vector  consisting  of  the  functions  a^,  32,  .  .,  a^.   This 
is  known  as  the  state  space  representation  of  the  system. 

PHYSICAL  CONSTRAINTS 
The  physical  constraints  on  the  state  variables  and  the  control 
inputs  have  to  be  specified  to  limit  the  state  space  to  the  required 
region.   Generally,   the  initial  and  final  states,   the  range  of 
controls  and  states,  and  time  may  be  specified. 
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THE  PERFORMANCE  MEASURE 
In  general,  the  performance  of  a  system  is  evaluated  by  a  measure 

of  the  form  [25] 

t^ 
J  -  h(x(tf),  tj)  +  ;^.  g(x(t),  u(t),  t)  dt,         (2.5) 

where  t^  and  t^  are  the  Initial  and  final  times,  h  and  g  are  scalar 
functions.   The  final  time  t^  may  be  specified  or  free. 

The  performance  measure  is  unique  for  each  set  of  control  and 
state  variable  values.  The  'h'  function  is  generally  used  to  specify 
constraints  to  be  satisfied  at  the  final  time.  The  'g'  function  is 
used  to  specify  the  time  varying  and  control  dependent  part  of  the 
performance  index.  The  performance  index  of  the  minimum  time  control 
problem,  without  constraints,  can  be  written  as 

J-L^ldt.  (2.6) 

THE  OPTIMAL  CONTROL  PROBLEM 
The  optimal  control  problem  is  to  find  the  admissible  control 

u*(t)  which  causes  the  system  described  by  the  set  of  first  order 

differential  equations  (2.4)  to  follow  an  admissible  traiectory  x  (t) 
that  minimizes   (or  maximizes)   the  performance  measure  defined  by 

(2.5).   The  control  u  (t)  is  called  an  optimal  control  and  x  (t)  is 
called  an  optimal  trajectory. 
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A  history  of  control  Input  values  which  satisfies  the  control 
constraints  during  the  entire  time  interval  from  the  initial  to  the 
final  time  [t^,  t^]  is  called  an  admissible  control. 

A  history  of  state  values  which  satisfies  the  state  variable 
constraints  during  the  entire  time  interval  from  the  initial  to  the 
final  time  [tp,  t^]  is  called  an  admissible  trajectory. 

CALCULUS  OF  VARIATIONS 
The  calculus  of  variations  is  a  branch  of  mathematics  that  is 
extremely  useful  in  solving  optimization  problems.  The  optimal 
control  problem  is  to  determine  the  control,  which  is  a  function  of 
time,  that  minimizes  a  performance  measure,  which  is  a  function  of 
functions  or  better  known  as  a  functional. 

FUNCTIONALS 

A  functional  is  best  described  as  a  function  of  a  function  or 
functions.  By  definition  [25],  "A  functional  J  is  a  rule  of 
correspondence  that  assigns  to  each  function  x  in  a  certain  class  Q  a 
unique  real  number.  Q  is  called  the  domain  of  the  functional,  and  the 
set  of  real  numbers  associated  with  the  functions  in  0  is  called  the 
range  of  the  functional . " 
For  example,  if 

yi"  fi^^r  ^^2^  ^^-^^ 

and  72  '  ^l^^V   "^2^  ^^"^^ 
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where  x,,   X2  are  independent  variables  and  f^,      ^2     ^^®   scalar 

functions,  then  the  quantity 

J  -  gCyp  72)  ^2.9) 

is  a  functional. 

VARIATION  OF  A  FUNCTIONAL 

The  variation  of  a  functional  is  analogous  to  the  differential  of 
a  function  in  that  it  helps  in  the  determination  of  the  extreme  values 
(maximum  or  minimum)  of  the  functional.  The  differential,  df,  of  a 
function,  f,  of  variables,  x^^,  X2 x^,  is  given  by 

df  -f£dx,  +!!dx,  +  .  .  .  +!£dx^-  (2.10) 

ax^     3X2  ax^ 

Similarly,  the  variation,  5J ,  of  a  functional,  J,  of  functions  y^,  y2, 
,  .  .  ,  y  is  given  by  the  relation 

6J-''  Sy.^'lsy^^  .   .   .  ^'-Sy^-  (2.11) 

3yi      3y2  dy^ 

FUNDAMENTAT.  THEOREM  OF  CALCUI-TTS  OF  VARIATIONS 

The  fundamental  theorem  states  that  the  variation  should  be  zero 
on  an  extremal  curve,  provided  there  are  no  bounds  on  the  curves. 
Mathematically  speaking,  if  x  is  a  vector  function  of  t  in  the  class 
fi,  and  J(x)  is  a  differentiable  function  of  x,  then 

«J(x*,  6x)-0  (2.12) 
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for  all  admissible  Sx.  An  admissible  Sx  is  one  which  satisfies  the 
condition  (X  +  5x)  e  n.  If  0  is  a  class  of  continuous  functions,  then 
X  and  5x  must  both  be  continuous. 

CONSTRAINED  MINIMIZATION  OF  FUNCTIONALS 
So  far,  we  have  discussed  functionals  involving  the  state  vector, 
x(t),  with  the  assumption  that  they  are  independent.  However,  this  is 
not  the  case  in  control  problems  where  the  state  trajectory  is 
determined  by  the  state  equations  (2.4).  This  involves  (n+m) 
functions,  x(t)  and  u(t).  of  which  only  the  m  controls  are 
independent.  The  state  vector  is  dependent  on  the  controls. 
Constrained  functionals  are  generally  minimized  by  the  Lagrangian 
multiplier  method. 

THE  LAGRANGIAN  METHOD  FOR  CONSTRAINED  MINIMIZATION 
The  Lagrangian  method  is  used  to  determine  the  minimization  of  a 
functional,  J,  of  the  form 

J(x,  u)  -  //  g(x(t),  u(t),  t)  dt,  (2.13) 

where  x(t)  is  the  state  vector  of  the  order  n,  and  u(t)  is  the  control 

vector  of  the  order  m,  subject  to  constraints, 

c(x(t),  u(t),  t)  -  0,  (2.14) 

where  c  is  a  vector  of  equality  constraints  of  order  n.   The  method 

consists  of  forming  an  augmented  functional 

t^ 
j(x(t),  u(t),  A(t),  t)  -  LMg(x(t),  u(t),  t) 

+  A'^(t)  c(x(t),  u(t),  t)]dt,     (2.15) 
21 


where  A.(t),  i-1.  2 n,  are  called  Lagrangian  multipliers .   The 

procedure  now  allows  us  to  find  the  function  u(t)  which  extremizes 
equation  (2.15)  and  satisfies  the  constraints  simultaneously.  Now,  if 
the  constraints  are  of  the  form  of  equation  (2.4).  then  the  augmented 
functional  becomes 

j(x(t),  u(t).  A(t),  t)  -;/{g(x(t),  u(t),  t)  +A'^(t) 

[x(t)  -  f(x(t),  u(t),  t)])  dt.   (2.16) 

When  the  constraints  are  satisfied  the  augmented  functional,  J,  equals 
the  unconstrained  functional,  J,  for  all  values  of  X(t)  and  time.  It 
should  be  noted  here  that  the  representation  of  the  state  equations  in 
equation  (2.16)  is  not  conventional.   If  we  define 

T 
g(x(t).  x(t),  u(t).  A(t),  t)  -  g(x(t),  u(t),  t)+A  (t) 

[x(t)  -  f(x(t),  u(t),  t)]       (2.17) 
then, 

J(x(t).  x(t).  u(t),  A(t).  t)-;/i(x(t).  x(t).  u(t).  A(t),  t)dt.  (2.18) 

0 

The  functional,  J,  can  be  minimized  by  setting  the  variation  5J 
to  zero.  The  variation,  SJ ,  after  simplifying  [Appendix  1]  is  given 
by 
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fij  -  fi  (x(t).  i(t).  u(t),  A(t).  t) 


t-t. 


dx(t^) 


-  [^     (x(t),  i(t),  u(t),  A(t).  t)  x(t)] 

ax 


dt. 


t-t. 


+  g  (x(t),  x(t),  u(t),  A(t),  t) 


t-tf  ''"f 


+  //{[!!  (S(t).  x(t),  u(t).  A(t),  t) 

^0  ax 


^  (!i  (x(t),  x(t),  u(t),  A(t),  t))]  5x(t) 

^^  ax 


+^  (x(t),  x(t),  u(t),  A(t).  t)  5u(t) 

au 


+  f!  (x(t),  x(t),  u(t).  A(t),  t)  5A(t))  dt.     (2.19) 

aA 

For  an  extremal  curve 

5J(x*(t),  x*(t),  u*(t).  A*(t).  t)  -  0.  (2.20) 

where  *  signifies  the  optimal  values.   For  a  minimal  time  control 
problem,  as  stated  before  by  equation  (2.6), 

g(2c(t),  u(t).  t)  -  1.  (2.21) 

Therefore , 

|(x(t).  x(t),  u(t),  A(t),  t)-l  +A'^(t)  [i(t)  -  f(x(t),  u(t),t)].  (2.22) 
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Substituting  equation  (2.22)  into  equation  (2.19)  produces 
5J  =  A''^(t^)  dx(t^)  +  [1  +  A'^(t)  (x(t) 


f(x(t),  u(t),  t))  -  A'^(t)  x(t)] 


t-t^     f 


+  A  {[-a'^CO  fz  (X(t).  M(t),  t)  -  i^(t)]  6s(t) 

^0      ax 

T 

+  [-A'^(t)  fz  (2^(t).  u(t).  t)]  5u(t) 
dn 

+[x(t)  -  f(x(t),  u(t),  t)]  5A(t))  dt.  (2.23) 

Since  6J  -  0  ,  equation  (2.23)  gives  us  the  necessary  conditions  for 
the  optimal  control,  which  are, 

S*(t)  -  f(x*(t),  u*(t),  t),  (2.24) 

i*(t)  -  -  f^^x*(t).  u*(t).  t)  A*(t),  (2.25) 

ax 

A*\t)  f£  (x*(t).  u*(t))  -0.  (2.26) 

and  1  -  A*'^(V  x*(tf)  -  0.  (2.27) 

The  n  equations  (2.24)  are  the  state  equations.  The  n  equations 
(2.25)  are  called  the  co-state  equations  or  the  multiplier  equations. 
The  m  equations  (2.26)  are  the  optimality  conditions  which  specify  the 
control  in  terms  of  the  state  and  costate  variables.  The  equation 
(2.27)   constitutes  the  boundary  condition  on  the  multipliers  or  the 
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transversality  equation.   The  transversality  condition  provides  one 
equation  for  the  unknown  final  time. 

The  above  conditions  are  not  sufficient  to  solve  the  time  optimal 
control  problem.  The  above  equations  would  lead  to  an  optimal 
solution  that  would  approach  zero  as  the  controls  move  towards 
infinity.  This  would  require  that  the  control  be  constrained.  In 
most  practical  systems,  however,  bounds  on  the  control  exist  in  the 
form  of  maximxam  force  or  torque  that  can  be  applied  by  the  actuator. 
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CHAPTER  3 

THE  CONTINUUM  APPROACH 

In  Chapter  2  the  general  time  optimal  control  problem  was 
formulated.  The  system  state  equations,  the  multiplier  equations,  the 
transversality  equation,  and  the  optimality  conditions  were  derived 
for  a  general  system.  However,  bounds  on  the  control  were  not 
specified,  whereas  physical  systems  have  bounded  controls.  In  this 
chapter  the  minimum  time  problem  with  bounded  control  will  be 
introduced  and  the  continuvun  approach  will  be  developed  for  the 
analysis  of  one  such  problem  in  the  phase  space  domain.  The 
solution  of  the  system  of  equations  derived  by  the  continuum  approach 
will  be  presented  in  Chapter  4. 

MINIMUM  TIME  MOTION  WITH  BOUNDED  CONTROL 
As  an  illustration  of  the  method  by  which  control  bounds  are 

treated  consider  an  n    order  linear  system  in  the  phase  variable 
canonical  form 

x(t)  -  A  x(t)  +  b  u(t)  (3.1) 

where  x(t)  is  the  state  vector  of  order  n,   u(t)  is  the  scalar  control 
input,  and  A  is  an  n  x  n  matrix  of  the  form 
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A  - 


~0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 
0 

1^1 

-^2 

■^3 

-^ 

-^5      •         • 

-a 

n 


(3.2) 


and  b  is  an  n  X  1  vector  given  by 


b  - 


0 
0 
0 


(3.3) 


0 
1 

The  control  magnitude  satisfies  the  constraints 

-1  <  u  <  1.  (3.4) 

The  initial  and  final  conditions,  x(0)  and  x(tj) ,  are  specified.   The 

system  described  by  equation  (3.1)  leads  to  an  augmented  functional, 
J ,  of  the  form 

J  =  Jq^  {1  +  A^(x  -  A  X  -  b  u)  +  t^(u   -    1)  +  n'(-u-l))    dt   (3.5) 

where  A  is  a  vector  of  Lagrangian  multipliers  of  order  n  and  /i  and  fj, 
are  non-negative  scalar  Lagrangian  multipliers  used  to  append  the 
control  magnitude  constraints  to  the  functional.   It  can  be  shown  [26] 

that  (i^     is  greater  than  zero  when  u  is  on  the  +1  boundary  and  zero 
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otherwise.  A  similar  statement  can  be  made  for  n'   with  respect  to  the 
-1  boundary. 

The  functional,  J,  can  be  minimized  by  setting  the  variation,  5 J , 
to  zero.   The  variation,  5 J ,  is  given  by 


5J  -  [1  +  a"^  (x  -  A  X  -  b  u)  +  /i"*"  (u-1)  +/i'(-u-l)] 


tf  ^^f 


^4 

^0 


+  fJ   (A^  fix  -  a"^  A  5x  +  (/  -  p'-  A  b)  5u 


+  (x  -  A  X  -  b  u)  6A  +  (u  -  1)  5^"^  +  (-U  -  1)  6fi')dt.    (3.6) 
The  second,   third,   and  fourth  terms  of  the  dt^  terms  go  to  zero  for 


all  values  of  t.   Integrating  by  parts,  and  using  5x(t^)  -  -  x  dt^ 


[26] ,  we  get 


6J  -  [1  -  A  x] 


+  SJ   d^   6s  -  a''^  A  5x  +  (m"^  -   t^'-  ^^   b)  Su 
tj    u 


+  (X  -  A  X  -  b  u)  5A  +  (u  -  1)  5;i  +  (-u  -  1)  Sti'}dt.  (3.7) 
If  we  choose  A  such  that  the  coefficients  of  6x  go  to  zero  and  since 
the  other  variations  in  equation  (3.7)  are  arbitrary,  we  get 


T  ' 
1  -  A^  X 


t-t. 


T 
A  A, 


0, 


/  -  /i"  -  A^  b  -  0, 


(3.8) 

(3.9) 
(3.10) 
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and                      i  -  A  X  +  b  u.  (3.11) 

The  last  two  terms  under  the  integral  sign  in  equation  (3.7)  indicate 
that  the  control  u  must  have  a  magnitude  of 

|u|  -  1.  (3-12) 

Equation  (3.10)  shows  that 

and  it   is   seen  that  the  variations  Sm"^  and  6^"  are  not  independent. 
Combining  equation   (3.13)  with  the  conditions  on  the  multipliers,  ^i^ 

and  /i  ,  we  get 

u(t)  -  sgn  (A^(t)).  (3.14) 

which  is  also  demostrated  in  optimal  control  texts  [25] -[26]. 

In  this  work  a  double  integrator  problem  has  been  considered. 
The  double  integrator  may  be  defined  as  one  in  which  the  control 
specifies  the  acceleration  of  the  system.  The  position  may  then  be 
calculated  by  integrating  the  control  twice.  The  double  integrator 
problem  may  be  defined  by  an  equation  of  the  form  (3.1),  where 

^  -  [2  J]  "'i^' 

and 

b-[;].  (3-16) 

Equations  (3.8),  (3.9),  (3.11),  (3.12)  and  (3.14)  give  us 

1  -  A^(t)  X2(t)  -  A2(t)  u(t)  -  0,  (3.17) 
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X^(t)-0,  (3.18) 

A2(t)  -  -  A^(t).  (3.19) 

x^(t)  -X2(t).  (3.20) 

X2(t)  -  u(t).  (3.21) 

|u(t)|  -  1.  (3.22) 

and  u(t)  -  sgn  (A2(t)).  (3.23) 

Equations  [3.8],  [3.9],  [3.11],  [3.12],  and  [3 . 14]  along  with  the 
initial  and  final  conditions  on  the  state  variables  provide  us  with  a 
problem  similar  to  the  two  point  boundary  value  problem  [TPBVP]  that 
can  be  solved  by  numerical  techniques.  However  we  lack  boundary 
conditions  on  the  multiplier  values. 

As  stated  in  Chapter  1  most  workers  used  some  sort  of  a  shooting 
technique  to  solve  the  problem.  However  most  of  these  techniques  are 
time  consuming  and  are  unfeasible  for  the  design  of  an  online 
controller.  Shetty  [19]  has  shown  that  conventional  numerical 
techniques  in  the  time  domain  are  unstable  unless  the  initial 
estimates  are  close  to  the  actual  solution.  The  various  novel  methods 
discussed  in  Chapter  1  either  require  a  linear  system  or  do  not  solve 
the  problem  completely.  The  switching  curve,  which  is  a  map  of  the 
sign  of  the  control  at  any  point  in  the  state  space,  is  the  fastest 
technique  and  still  the  best  control  method  but  is  not  possible  for 
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higher  order  systems  because  it  is  hard  to  visualize  and  only  of 
limited  use  for  nonlinear  systems. 

THE  CONTINUUM  APPROACH 
Luh  and  Shafran  [24]  developed  a  method  to  determine  the  control 
from  a  least-squares-fit  of  the  distribution  of  the  minimum  time 
isochrones  in  the  phase  space.  The  isochrones  are  surfaces  of 
constant  final  time.  This  type  of  an  analysis  provides  a  convenient 
method  to  determine  the  control  for  an  infinite  number  of  different 
problems  in  an  efficient  manner.  Linear  systems  require  a  single 
distribution  of  isochrones,  however,  nonlinear  systems  require  one 
isochrone  distribution  for  each  exclusive  destination.  Whereas  Luh 
and  Shafran  determined  the  distribution  from  many  different  solutions 
of  a  linear  system,  the  goal  in  this  work  is  to  determine  if  this 
isochrone  distribution  could  be  determined  from  some  balance  relation 
involving  the  gradients  of  the  final  time. 

There  are  several  points  which  promote  this  idea  and  provide  an 
incentive  for  the  work.   These  are  [27] 

1.  Any  exclusive  distribution  is  a  function  of  the  state 
variables  exclusively,  thus  eliminating  time  as  the 
independent  parameter. 

2.  Since  the  state  variables  are  the  independent  parameters,  the 
steps  in  finding  the  isochrones  are  the  same  for  both  linear 
and  nonlinear  systems. 
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3.  If  the  isochrone  information  can  be  determined  quickly  then  it 
may  be  possible  to  develop  a  controller  based  on  these 
principles. 
A  continuum  approach  in  the  phase  plane  is  presented  for  a  double 

integrator  problem. 

GOVERNING  EQUATIONS 

The  necessary  and  sufficient  conditions  for  a  minimum  time 
control  problem  are  given  by  equations  (3.8),  (3.9),  (3.11),  (3.12) 
and   (3.13).   The  co-state  vector  A  at  the  initial  time  is  related  to 

the  final  time  by  [26] , 

A(t^)  -  -  Vtj(s(tj^)).  (3.24) 

where  V  is  the  vector  gradient  operator  given  by 

rT  -  (  3   a   _  _  a   j^  (3.25) 


r  -  ( 


ax^  3x2      dx^ 


for  an  n^^  order  system.  Equation  (3.25)  holds  over  all  the 
permissible  regions  in  the  phase  space  where  the  gradient  exists. 
Equation  (3.24)  will  be  derived  for  a  general  system  defined  by 
equation  (2.4). 

The  augmented  functional,  J,  for  a  minimum  time  control  problem 
driving  a  system  to  the  origin  can  be  written  as 

J  -  SJ   [1  +  k^   (^  ■   £(x,ii))]  dt  -  J(x(t.)  -  tj(x(tj^).(3.25.1) 

This  is  a  continuous  function  of  s(tj^).  Integrating  equation  (3.25.1) 
by  parts  produces 
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-j(x(tj^)  -Ax 


^f  +  /^  [1  -  a"^  X  -  a"^  f(x,u)]  dt.   (3.25.2) 


Applying  the  gradient  operator  then  produces 


V  J  -  7t^(x(t.)  -  -A^Ct.)  +  L^  Y  [1  -  A^  X  -  A^  f(x.u)]  dt. (3.25.3) 


Since  V  consists  of  partials  with  respect  to  the  components  of  x(t^) 
then  the  integrand  in  equation  (3.25.3)  vanishes  for  all  t  not  equal 
to  t..   For  the  case  when  t  is  equal  to  t^^  we  have 


V[l  -  A^  X  -  a"^  f(x,u)] 


A*^  -  A*^  V  f(x,u) 


(3.25.4) 


Rearranging  equation  (3.25.4)  produces 

A  +  [V  f(x,u)]'^  A  -  0.  (3.25.5) 

which  shows  that  the  integrand  of  equation  (3.25.5)  vanishes  for  all  t 
in  the  interval  (0.  t j) .  Hence  we  get  equation  (3.24)  from  equation 
(3.25.3).  It  should  be  noted  that  due  to  the  choice  of  signs  of 
multipliers  in  equation  (3.25),   the  vector  A(t^)  points  in  the 

direction  of  the  greatest  decreasing  final  time. 

The  transversal ity  condition  equation  (3.8)  is  usually  applied 
only  at  the  final  time,  but  it  is  true  everywhere  in  the  phase  space 
where  the  gradient  exists.  The  transversal ity  condition  evaluated  at 
the  initial  state  is 


1  +  S  (t^)  Vtj(x(tj^))  -  0. 


(3.26) 
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Equation  (3.26)  is  the  projection  of  the  gradient  of  the  final  time  on 
the  time  derivative  of  the  state  vector.  Equation  (3.26)  is 
insufficient  to  uniquely  specify  the  components  of  the  gradient  of  t^ 

for  systems  of  order  two  or  greater. 

NECESSARY  CONDITIONS 

In  order  to  develop  a  method  to  determine  the  final  time  field, 
t  (x(t.)),   it  is  necessary  to  develop  a  means  of  uniquely  specifying 

the  gradient  of  the  final  time. 

The  cost  function  of  the  minimum  time  problem  is  given  by 

t, 

-f-J 

which  can  be  written  in  the  augmented  form  as 


-Jo^ldt,  (3.27) 


^f  -  V  (1  -  ^  2C  +t   X  )dt.  (3.28) 

It  should  be  noted  that  the  final  time  value  is  the  same  as  the 
augmented  functional.  The  first  two  terms  of  the  integral  constitute 
the  transversality  equation.   Equation  (3.28)  then  reduces  to 

tf  -Jq^  A^xdt.  (3.29) 

As  stated  earlier  in  this  chapter,  the  double  integrator  will  be  used 
to  illustrate  the  problem.  Integrating  equation  (3.29)  by  parts  for 
the  double  integrator  and  using  state  and  costate  equations,  produces 


T 
tf  -  A  X 


o'  -  Jo^  A^xdt, 
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or 


tf  -  A  X 


0^  -  Jo^  (Xi  '\   +  X2  X^)  dt. 


The  first  term  under  the  integral  is  zero,  and  the  second  term  under 
the  integral  can  be  rewritten  using  the  state  and  costate  equations. 
This  gives 


T 
tj  -  A  X 


o'  ^  J"o'  ^1  ^1  ''^- 


Integrating  again,  we  have 

t/ 


tf  -  A  X 


0  -^  ^1  ^1 


j/  X^   x^  dt. 


The  term  under  the  integral  is  zero  by  equation  (3.18).   Therefore,  we 

have 

^f 
tj  -  (2  x^  X^  +   X2  X^)      0  •  ^^-^^^ 

If  the  destination  is  the  origin  of  the  phase  plane,  equation  (3.30) 

becomes 

t^  -  -  2  Xj^(O)  A^(0)  -  X2(0)  A2(0), 

which  gives 

tj  -  [2  x^(0)   X2(0)]  7t^(x(0)).  (3.31) 

Equation  (3.31)   is  the  additional,  necessary  condition  to  uniquely 
specify  the  gradient. 

As  stated  before,  the  minimxim  time  problem  is  solved  in  the  time 
domain  approach  by  the  integration  of  the  functional  along  the  phase 
trajectory  and  the  functional  value  is  the  difference  between  the 
final  time  and  the  initial  time,  i.e.  the  endpoints  of  the  path  of 
integration  are  the  initial  and  final  points.   Equation  (3.30)  shows 
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that  the   final  time  can  be  described  as  a  function  of  the  local 
variables  at  each  end  point.   Now,  equation  (3.29)  can  be  written  as 

tf  -  Jo^  ^  k   dt, 


or, 


or. 


t.-j: 


X(tj)     ^ 


x(0) 


A  dx, 


x(t.)   „ 

-J^(0)  ^%^' 


(3.32) 


which  shows  that  t^  is  a  potential  function.   The  path  of  integration 
is  arbitrary  and  the  value  depends  only  on  the  endpoints  of  the  path 

of  integration. 

It   is   instructive  to  solve  equations   (3.26),   and  (3.30) 
simultaneously  for  the  components  of  the  gradient  of  t^.   The  solution 


IS 


^^f 


(2  x^  u  -  X2) 


u 

-X, 


2  X, 


^f 
-1 


(3.33) 


"2    "     1 

where  it  can  be  seen  that  the  coefficient  matrix  is  singular  on  the 
switching  curve  [26].  Equations  (3.26)  and  (3.30)  are  used  to 
determine  the  final  time  field.  These  equations  are  solved  by  a 
finite  element  method  [28].  The  details  of  the  finite  element  method 
are  given  in  Chapter  4. 
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CHAPTER  4 

FTNTTR  ELEMENT  MODEL 
INTRODUCTION 
Engineering  problems  are  invariably  nonlinear  in  nature  and 
closed  form  solutions  are  not  available.  This  necessitates  the  use  of 
some  form  of  numerical  technique.  Two  methods  which  are  commonly  used 
are  the  finite  difference  and  the  finite  element  methods.  The  finite 
element  method  uses  piecewise  continuous  approximations  of  the  region 
under  consideration.  The  finite  element  method  is  more  versatile  and 
easier  to  use  with  problems  that  have  irregular  geometry  or  unusual 
boundary  conditions.  The  finite  element  method  was  used  in  this  work 
because  of  the  success  of  the  method  in  treating  other  problems  based 
upon  a  continuum  formulation. 

The  finite  element  method  discretizes  the  solution  region  of  a 
continuum  problem  into  a  finite  number  of  elements.  The  unknown 
solutions  within  each  element  are  then  expressed  in  terms  of 
approximating  functions.  The  problem  is  reduced  to  one  with  a  finite 
number  of  unknowns,  the  nodal  values,  and  solved. 

The  advantage  of  such  a  formulation  is  that  it  reduces  the 
complex  problem  to  a  greatly  simplified  problem  at  an  elemental  level. 
The  method  provides  a  variety  of  different  ways  in  which  the  elemental 
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problem  can  be  formulated.   The  least-squares  method  has  been  used  in 
this  work  to  produce  the  element  equations. 

The  finite  element  method  has  five  basic  steps.   These  are, 

1.  Discretization  of  the  continuum. 

2.  Specification  of  the  interpolation  function. 

3.  Development  of  the  element  equations. 

4.  Assembly  of  the  element  equations  to  get  the  system  of 

equations . 

5.  Application  of  the  boundary  conditions  and  obtaining  the 
solution  of  the  system  of  equations. 

In  this  Chapter  the  finite  element  model  of  the  double  integrator 
problem,  for  which  the  continuum  relations  were  derived  in  Chapter  3, 
will  be  presented.  The  presentation  of  the  results  and  the  pertinent 
discussions  and  conclusions  will  be  deferred  until  Chapter  5. 

FINITE  ELEMENT  ANALYSIS 
As  stated  before,  the  finite  element  method  basically  consists  of 
five  steps.  In  this  work  the  analysis  was  done  with  the  help  of 
FORTRAN  programs  executed  on  both  an  IBM  370  VM/CMS  computer  and  an 
IBM  PC  microcomputer.  The  program  provides  for  changes  in  the  size  of 
the  region  under  consideration  and  the  size  of  the  grid,  but  it  allows 
only  a  rectangular  grid  with  triangular  elements.  The  least  squares 
technique  was  used  in  this  analysis  so  as  to  provide  symmetrical 
element  matrices  for  economical  storage  and  faster  solution  time.  The 
analysis  was  performed  using  single  precision  arithmetic  and  is 
suitable  for  implementation  on  microprocessors.   The  steps  in  the 
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formulation   and  analysis  of  the  finite  element  model  will  be 
presented.   The  programs  are  listed  in  Appendix  2. 

DTSCRETIZATTON  OF  THF.  CONTINUUM 

The  discretization  of  the  continuum  has  been  done  with  two 
FORTRAN  subroutines,  namely,  NODE  and  GRIDLT.   The  domain  selected  is 
a  square  region  with  the  origin  at  the  center.   The  subroutine  NODE 
places  equispaced  nodes  in  two  directions  and  numbers  them.   The  nodes 
are  numbered  horizontally  if  the  number  of  nodes  in  the  y-direction  is 
greater  than  or  equal  to  the  number  of  nodes  in  the  x-direction  or 
vertically  downwards  otherwise   (Figure   (4.1)).   This  reduces  the 
memory  required  when  the  element  equations  are  stored  in  banded  matrix 
form.    It  should  be  noted  that  there  is  no  overall  advantage  to  any 
method  when  the  number  of  nodes  is  the  same  in  both  directions.   The 
nodes  divide  the  region  into  a  regular  and  uniform  grid.   The 
subroutine  GRIDLT  divides  the  grid  into  triangular  elements  of  the 
same  size  and  forms  a  table,  called  the  node  table,  which  contains  the 
node  numbers  of  each  element.   The  node  table  has  its  elements 
specified  in  the  counterclockwise  direction.   This  is  shown  in  Figure 
(4.2).   The  triangular  element  was  chosen  for  ease  of  application. 

INTERPOTATION  FUNCTION 

The  next  step  is  to  choose  the  interpolation  function.   The  nodal 

values  of  the  final  time,   t.,  are  t.  .  t.  ,  and  t^  and  the  nodal 

^       ^i    j       k 

coordinates  are   (x,  ,   x,  ) ,   (x,  ,   x,  ) ,  and  (x^^  •  ^j  >  •   A  li'^ear 

^i    ^i      j     j         k    k 
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I  2  3 


NX  NX+l 


CNX+0*l 


2CNX*U*l 


NY(NX*l)+l 


2(NX*U 


3<NX*U 


<NY*l)0«*l) 


NX   °   Nunbvr   of   Rvctangutar   Regions  n   th«   x   axis 
NY   =   Nunbvr   of   Rvctangutar      Regions   »   th*   y  axis 

Figure    4.1a;   Nunbering    of    Nodes 
Case    M   NY    1   NX 


NY*1 


CNY*l>*l    aCNY->0»l 


MX<NY*U*1 


2tNY*l)       3CNY*0  (NX*IXNY*U 

NX   =   Nunbrr   of   rectangular  regions  in   the   x   axis 
NY   -   nurtoer   of   rectangular   regions  •>   the   y   axis 


Figure    4.ro:    Nunbering    of    Nodes 
Case   \u    NY   <   NX 
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Figure    4.2a:    Basic    Element    Shape 


FiQ'^re    4.2b:    Nodal    Values    of    Final    Tine 
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interpolation  formula  [28]  of  the  form 

t^  -  a,  +  02  Xj^  +  a^   X2 

is  chosen.   The  nodal  conditions  give 

^f.  "  "1  "*■  "2  ''l.  ^  "3  "^1.' 


t^  -  a,  +  a„  X,   +  a.,  x„  , 


tc     -  on  +  "o  X,   +  a,  x«  , 
^k    ^    "^   k    "^  ^k 


(4.1) 


(4.2) 


which  yields 


"1  - 


2  A 


[(Xt  x„  -  x^  X,  )  t.  +  (X.,  Xj  -  x^  X2  )  t^ 
Ij  \       \   2j    fi    Ik  2i   li  ^k   ^j 


+  (x,  X„  -  X,  x„  )  t^  ] , 
^i  ^j   ^j  ^i   ^k 


a  -  _]_  [(x,  -  X2  )  t.  +  (X2  -  X2  )  t^  +  (X2  -  X2  )  t.  . 

2   2  A    ^j  ^k   ^i  ^k   ^i   ""j     ^1    J    k 

"3  -  J_  [(^1  -  x^  )  tf  +  (x^  -  X,  )  tf  +  (x^  -  x^  )  t.  , 

■^   2  A    ^k    J   ^i  ^i   "-k   "-J     J    1    k 


where  the  determinant 


1      X,      X2 
i       i 


^      ^1.     ^2, 
J       J 


1      X,       x„ 
^k      "^k 


2  A, 


(4.3) 


and  A  is  the  area  of  the  triangle. 
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Substituting  for  a^,  a^,    a^   in  equation  (4.1)  and  rearranging,  we 

get 

t.  -  N.  t.  +  N  t.  +  N,  t.  ,  (4.4) 

f    1   f^    J   tj     K  Ij^ 


where 


N.  -  ^  [a.  +  b.  X,  +  c.  x,  ] ,  (^-5) 

^   2  A 

N.  -  J_  [a.  +  b.  X,  +  c.  x,  ],  (4-6) 

J    2  A   -'    ^  ^ 

2  A 


and 


N. 

1 


-i  -  -1.  %-  X  '2j'  ^  -  "2j-  %•  ^i  -  V  "^j' 

N.    and  N,  are  called  the  shape  functions.   It  may  be  noted  that 
J'       k 

each  shape  function  is  one  at  its  node  and  zero  at  the  others,  the  sum 
of  the  three  shape  functions  at  any  point  is  one,  and  the  shape 
function  varies  linearly  between  its  node  and  the  other  two  nodes  but 
is  zero  on  the  side  opposite  its  node.  The  selection  of  this  type  of 
an  interpolation  function  gives  us  the  final  time,  t^,  as  a  linear 

function  of  x,  and  x^  and  also  gives  constant  first  derivatives  in  all 

elements.   The  first  derivatives  of  t^  are  given  by 
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and 


^^   =  ^   [b.  t.  +  b,  t.  +  b,  t.  ],  (4.8) 

^  -  1   fc.  t-  +  c.  t.  +  c,  t.  ].  (4.9) 


2  A       i    -^   J        k 


The  element  nodal  coordinates  can  be  transformed  into  a 
dimensionless  local  coordinate  system.  Such  a  transformation  is 
useful  only  in  that  it  provides  us  with  a  integral  relation  involving 
the  shape  functions.  The  details  of  the  transformation  are  given  in 
reference  [28]  and  it  gives  an  integral  of  the  form 

f  nPn^N^  dA  -    P'  '^'  ^'    2  A.  (4.10) 

"'A  1  J   k         (p+q+r+2) ! 

This   integral   is  useful   in  the  determination  of  the  element 

properties. 

ELEMENT  PROPERTIES 

The  element  properties  have  been  derived  with  the  least- squares 
technique.  Using  equations  (3.24)  and  (3.29),  we  seek  the  solution, 
t^(x-,,  X2)  ,  that  minimizes  the  integral 


I..  -  L  ( J_ 


fe   •'A 


(1  +  x^  Vtj)^  +  K  (tj-  [2x^   X2]  Vt^)^ 


+  ^i^(n   -  1)  +  /i'(-u  -  D)  dA,  (4.11) 

where  I^   is   the  performance  index  measuring  the  error  in  the 

solution  and  K  is  a  constant  of  unit   magnitude   that  insures 
consistency   of  dimensional  units.   Other  formulations  were  tried, 
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however,   equation   (4.11)   produced  the  results  which  fit  the 
requirements  and  by  which  the  solution  could  be  obtained. 

The  control  can  be  determined  to  minimize  the  performance  index, 
I   .    Setting  the  partial  derivative  of  I^^  with  respect  to  the 

control,  u,  to  zero,  we  get 

/i+  -  A2(l  -  X2  A^)  -  A^  ;  u  -  1,  (^-12) 

and  m"  -  -^2(1  -  X2  X^)  -  A^  ;  u  -  -1.  (4.13) 

Since  the  multipliers,  /i'^and  /i' ,  are  always  non-negative,  we  have 

A2(l  -  X2  A^)  >  A2  ;  u  -  1,  (^-1^) 

and  -A2(l  -  X2  A^)  >  A^  ;  u  -  -1.  (4.15) 

Equations  (4.14)  and  (4.15)  can  be  combined  to  give 

u  A2  (1  -  X2  A^)  >  0.  (4.16) 

Solving  equation  (4.16)  for  u  produces 

u  -  sgn  (A2  (1  -  X2  A^)).  (4.17) 

Equation  (4.17)   is  a  discrete  form  of  the  control.   The  control,  u, 
can  change   for  different  values  of  X2  within  the  element  and  the 

multipliers.   A,   and  A2,   are  constant  inside  each  element  but  vary 

between  elements.   This  gives  us  a  distinct  control  at  each  point  in 

the  phase  space  region  chosen.   The  control  is  chosen  such  that 
equation  (4.11)  is  always  reduced. 

Now,  expanding  the  first  term  of  equation  (4.11)  gives 
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^2 


X2  x^u 


x^u   u 


dA  Vtj,  (4.18) 


where  the  integrals  are  computed  over  the  area  of  the  element.  Since 
the  control  may  either  be  constant  or  vary  inside  each  element,  the 
integrals  involving  the  control,  u,  need  to  be  determined  separately. 
If  the  control,  u,  is  constant  within  the  element  then  the  integrals 
are  easily  evaluated.  However,  if  the  control  changes  sign  within  the 
element  then  the  integration  requires  another  procedure.  Assuming 
that  the  distribution  of  the  final  time,  t^,  is  known,  it  is  possible 

to  determine  whether  there  is  a  sign  change  within  the  element  by 
evaluating  the  control  equation  (4.17)  at  each  node.  If  there  is  a 
sign  change,  then  at  some  value  of  X2  the  argument  of  equation  (4.17) 

must  vanish.   This  value  of  x,  will  be  denoted  as  X2  and  is  given  by, 

s 

X,  -J_-__li__  (^-15) 

''s    A^   atj  /  dy.^ 

which  is  a  function  of  the  nodal  values  of  the  final  time.  Figure 
(4.3)  shows  the  two  possible  distribution  of  the  control,  u,  inside  an 
element  in  this  case.  It  should  be  noted  that  the  grid  was  completely 
made  up  of  these  two  types  of  elements  exclusively. 

There  are  two  integrals  in  equation  (4.18)  that  depend  upon  the 
control,  u.  With  reference  to  Figure  (4.3) 

;^  u  dA  -  u.  4^dA  -  u.  4.^^dA  -  u.(2A,  -  A).       (4.20) 
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J 


Figure    4.3a:    Division    of    Triangle    for    Integration    Purposes 
Case    i 


Figure    4.3to'    Division   of    Triangle    for    Integration   Purposes 
Case   ii 


47- 


The  area  A,  is  a  function  of  the  location  of  the  switching  line  given 
by  equation  (4.19)  and  is,  thus,  a  function  of  the  nodal  values  of  t^. 
The  area  A^  is  given  by 

A,  -  A  (X,   -  X.  )/(x„   -  Xj    )  (4.21) 

^       ^s  ^j    ^k    '^j 

It  should  also  be  noted  that  the  quantity  u^  in  equation  (4.20)  is  the 
nodal  control  value  that  has  a  unique  sign.  The  other  integral  in 
equation  (4.18)  that  depends  upon  the  control,  u,  is 

;^  X2  u  dA  -  u.  J^^X2  dA  -  u.   !^.^^^2   ^' 

4  X2  u  dA  -  2  u.  /^^X2  dA  -  u .  ;^  X2  dA.         (4.22) 

Now,  each  of  the  integrals  in  equation  (4.22)  is  the  first  moment  of 
area  of  the  triangles  about  the  x^  axis,  therefore,  equation  (4.22) 

becomes 

r,  x^  u  dA  -  u.[2  A,  (X,  +  2  x„  )  -  A  (x,  +  X2  +  X2  )]  /  3   (4.23) 
JA   2         J     i    2j      Zg        ^^        ^j    ^y. 

where  x„   is  the  X2  coordinate  of  the  node  having  the  control  u^  while 
J 

x^   and  X,   are  the  x„  coordinates  of  the  remaining  nodes.   The 
^i        \ 

remaining  integrals  in  equation  (4.18)  give 

r  x„  dA  -  A  (x„  +  X,  +  X,  )  /  3,  (4.24) 

Ja  2         2^       2j    2^ 

and  ;^u2  dA-4  dA-  A.  (4.25) 

Defining  the  vectors  and  the  matrix,  d,  as 
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Now 


N^  -  [  N^  Nj  Nj^  ] 


t^  -  [  t^  tf  tf  ]  , 


^     ^i  ^j  ^k 


-2 


[  X^   X,^   Xrt  J  , 

^i  ^j  ^k 


(4 

26) 

(4 

27) 

(4 

28) 

(4 

29) 

and 


d^. 


\       ^j   \ 
^i   ^j   '^k 


(4.30) 


we  have 


J^  x2  dA  -  ;^  N^  X2  N^  X2  dA, 


or 


;  x^  dA  -  ;^  4  N  N^  X2  dA, 


or 


J^xl   dA  -  4  Ja^  ^  ^^2' 


(4.31) 


;.  N  N*^  dA  -  ;, 


Ni 


N,N2 
N^N3 


N,N2 

^^2 
N2N3 


N1N3 
N2N3 

N? 


dA. 


Using  equation  (4.10)  we  have 


;^NN^ 


Defining  the  matrix  G  as 


dA  - 


"    2 

1 

I 

A 

- 

1 
1 

2 

1 

I 

12 

2 

2 

1 

1    " 

1 

2 

1 

t 

1 

1 

2 

(4.33) 


(4.33) 
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then  equation  (4.31)  gives 


2  ,,     A  .T 


;^  x^  dA  -  _2.  x^  G  x^. 


Now,  the  second  part  of  equation  (4.11)  gives 

2  X, 


(4.34) 


1  J. 


1  oT. 


Ja  -L  t^f  dA  .  _^  Y^t,  /, 


r  2  x^   Xj  1  dA  Vtj 


4  [  2  ^1  '^f   ^^2  ^f  ]  ^  ^^f 


or  I 


11 


-  r    ^  t^ 

fe   Ja  ^-  ''f 


t;  dA  +  Jl  v^t^  ; 

2 


-f  Ja 


4  X,    2  x^  X2 

2 
2  Xj^  X2     X2 


dA  Vt 


/a  [  2  ^1  '^f   ^^2  ^f  ] 


dA  Vtj. 


The  integral  in  the  first  term  of  equation  (4.35)  gives 
4  4_  t2  dA  -  _L  /^  n'%  N%  dA. 


-"f 


(4.35) 


or 


or 


/;,  i.  t|  dA  -  1,  4  4  H  H^  dA  t,, 


;.-14dA-i. 


'A  -2-  "f 


tf  G  tf . 


2    12 
The  integrals  in  the  second  term  of  equation  (4.35)  give 

J^  4  xj  dA  -  4  J^  N^^^  N^X]^  dA, 


or 


or 


and 


;^  4  x^  dA-  4  xj  ;^N  N^  dAx^. 
4  4  x2  dA  -  4  _^  x^  G  X,. 

4  2  x^  X2  dA  -  2  ;^  n\  N^X2  dA, 


(4.36) 


(4.37) 
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or 


or 


J^2  x^  X2  dA-2  x^J^N  N^  dAx2, 


4  2x,  X2dA-2_A_x{Gx2. 


The  integrals  in  the  third  term  of  equation  (4.35)  give 
;^  2  x^  tf  dA  -  2  J^  f  tf  e\   dA. 


(4.38) 


or 


or 


and 


;^  2  x^  tj  dA  -  2  4  ;^  N  N^  dA  x^. 


4  2  ^1  -f 


t^  dA  -  2   ^   tj  G  X,  , 
12 


J"a  ^^2  ^f  ^  ■  -^A  ^^   %  ^\   ^• 


(4.39) 


or 


or 


/, 


^  X2  tj  dA  -  tj  J"^  N  N  dA  X2, 


/a  X,  t^dA-_^  t^Gx2. 


12 


(4.40) 


Now,  from  equations  (4.8)  and  (4.9) 


and 


Vtj 


1^1 


d  tj, 


2  A 


t^.  dj. 


2  A 


(4.41) 
(4.42) 


If  we  define 


^1  "  4  ^""2   ""^  '^• 


X2U   u 


(4.43) 


(4.44) 


L 


4  X,    2  X,  X2 
2 


dA. 


(4.45) 
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and 


^4- 


12 


2  G  X,    G  ^2 


Then,  equations  (4.18)  and  (4.19)  can  be  combined  to  give 
I^  -  ^  +   ^   P.  /  t.  +  ^  _]_  4  ^  ^2  ^^  if 


(4.46) 


12 


^   4  A"^ 


^   £f  £4  ^"^  if 


2  A 


(4.47) 


an 


In  order  to  minimize  the  error  in  the  element  we  have  to  set  the 
partial  of  I^^  with  respect  to  the  nodal  values  t^  to  zero.   This  c 

be  done  by  making  use  of  the  vector  differential  properties. 

,T 
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fe 


atj 


-  0. 


which  gives 


1      T     1 
d  Pj  + 


2  A 


2  A 


'^   +   ^   d  £2  d'''  t^ 


+   ^   G  t.  +   •*•   d  P,  d'^  t .  +  _ 


4A' 

1    1 


4A 


^T  ,   aPo   ^T. 
tir  d    2   d  t^ 

atj 


1       T 

£4^  % 


^   d  pT  t^  -  0, 


(4.48) 


2  A   ""     ^    2  A 
The  second  term  in  equation  (4.48)  gives 


'h     d^t. 


2  A 


atj 


2  A 


1  d  t^  


ax. 


atj 


(4.49) 
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Evaluating   the   partial   differential  terms  in  equation  (4.49) 
separately,  we  have 


and 


3i. 

axo 


(Xo  -  X,  ) 


2,    2. 


ax. 


atj 


(  -  1/  a:^  )  (■ 


J  J 

1 
TX 


)  b' 


(4.50) 


or 


ax 


atj 


s  -  Cj^  b  , 


(4.51) 


where 


2  A  a: 


(4.52) 


Equation  (4.49)  can  be  rewritten  as 


'h     d^t, 

2  A   at^ 


2  A 


Cj^  b  P5  d  t^ 


(4.53) 


where 


^  3x2 


(4.54) 


The  sixth  term  in  equation  (4.48)  gives 


1    1 


2     2 


t^r  d    2  d  t^ 

atj 


Til 


2   4a2 


^T  ,  ap,  ,T. 

tf   d    I     d  t^ 

^   ax.. 


ax. 


at^ 


^JL_[.^£e^^%^i^-^r. 


2   4a2 


J^ 1_  c^  b  tj  d  Pg  d"^  tj,   (4.55) 


2   4a2 
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where 


P  -  ^^2 

^   axo 


or 


^6 


dv 


2(1.2) 


ap 
ax. 


ax. 


2(1.2)     0 


(4.56) 


The  partial  differentials  in  equation  (4.56)  can  be  evaluated  as 


^^2(1.2)   -  "i 


ax. 


2  A 


(X2  +  X2  ) 
(x,   -  x„  ) 


+  4  A, 


2  u 


or 


X 


(X2  +  X2  ) 
j     s 
(Xo   -  x«  ) 


+  2  A, 


(4.57) 


Equation  (4.48)  can  now  be  rewritten  as 


[EM(t.)]  t.  +  V(t.)  -  [0  0  0]'^  -  [RES]^ 


(4.58) 


where 


[EM(t^)]  - 


^   c^bP^/  +  _L    dP.  /+   ^ 


2  A 
1 

77 


-  -2  - 


4  A' 


12 


d  Po  d^-  _^  P4  d^ 
^      2  A 


1      T 
^   d  P, 
-  —4 


(4.59) 


2  A 


and 


[V(t  )]  -  _J_  d  P{  +  jL  _^: c^  b  t^  d  Pg  d^  tf   (4.60) 


2  A 


2   /  a2 
4  A 


This  gives  us  three  elemental  equations  for  the  three  nodal 
unknowns.   However,   the  nodal  values  are  also  associated  with  other 
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nodes  and  the  problem  has  to  be  solved  on  a  global  basis.   The  element 
equations  are  computed  in  the  FORTRAN  routine  LINTRI. 

ASSEMBLY  OF  THE  ELF-MF.NTAT.  EQUATIONS 

In  order  to  solve  for  the  unknown  nodal  values  of  the  final  time 
in  the  region  under  consideration.  It  is  necessary  to  combine  or 
assemble  the  element  matrix  equations  to  form  the  global  system  of 
equations.  The  basis  for  the  assembly  procedure  is  that  all  the 
elements  are  interconnected  at  the  nodes  with  adjacent  elements  and 
the  value  of  the  unknowns  are  the  same  for  each  element  sharing  that 
node.  The  assembly  procedure  is  carried  out  by  means  of  the  routine 
BUILD. 

SOLUTION  OF  THE  SYSTEM  OF  EQUATIONS 

In  order  to  solve  the  system  of  equations,  the  values  of  final 
time,    t^,   were  specified  as  boundary  conditions  on  the  outer 

boundaries.  This  was  performed  in  routine  BC.  The  solution  was  then 
obtained  iteratively  by  using  Newton's  method.  In  order  to  solve  the 
system  of  equations  the  Jacobian  has  to  be  calculated.  The  Jacobian, 
[J],  can  be  computed  from  equation  (4.58)  and  is  of  the  form 

[J]  -  [EM(tj)]  +  _f_  [EM(tj)]  %  +  _i_  [V(tj)].     (4.61) 

The  second  term  in  the  equation  (4.61)  gives 
'         [EM(t,)]t,  -   ^   bP^dTt._!fl^  +  _i_d  J^d^  t, 

-IT,  '       '     ^r^       '        'at,        ,  ^2        at. 
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A_bPU^£f^^  bT  +  ^_d  J^d^  t^c^b^ 
2  A  ^  AA^        ^  ^2     3x2 


^   b  P^  /  t.  _^  b^  +  _1 d  Pg  d^  tf  c^b^.  (4.62) 

TT-   -5        AA^        ^^2 


The  third  term  of  equation  (4.61)  gives 


a 


dt^ 


[V(tf)] 


1   azi  ^  1   1 


2  A   dt^ 


2   /  a2 
4  A 


2  b  t^  d  Pg  d^  c-,^ 


+  b  t^  d  ^^6  d^  t.  c,  +  b  t'^  d  P,  /  t,  ±1 

"  "f  "  "atT "  ^  ^       ^    ^  at. 


or 


1   d  P,  c,  b^  -.  _i_  ^_ 
-l-T--'     ^  2   ,^2 

+  b  t^  d  ^-6  d"^  t^,  c,    ^ 


2  b  t^  d  Pg  d'^  c^ 


ax. 


'-f   "1 


atj 


+  b  t''^  d  P,  d''^  t^,  ^1   b'^ 
-  -  -  -6  -   f  A  A, 


(4.63) 


The  third  term  in  equation  (4.63)  can  be  evaluated  as 
^    ^    b  J  d  i!6.  d^  t,  c,  ^ 


2   /  *2 
4  A 


ax. 


'f  "1 


atj 


L   ^    b  t^  d  P^  /  tj  cj  b''^ 


2   /  a2 
4  A 


(4.64) 


where 


'       ax^ 
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or 


0 

2 

^^2a.2) 

0   2 

2 

ax^ 

£7- 

^  ^(1.1) 

2s 

0   2 

0 

ax2 

L       ^ 

J 

The  non-zero  terms  in  equation  (4.65)  can  be  evaluated  as 


(4.65) 


3^^222  -  «  u. 


ax 


3   ^   (x^   -  X,  ) 


2  'k 

s 


(4.66) 


Equation  (4.63)  can  now  be  written  as 

^    [V(tj)]  -  \      d  £5  c^  b*^  +  J \ 


at 


2  A 


2    4a2 


2  b  t W  P,  d  c-i 


T       T     2   T 
+  b  t^  d  P^  d  %  Ci  b 


T         T      C      T 

+  b  t  d  P,  d  tr:   1^  b 

-  -     -     ^  *   A  A, 


(4.67) 


The  Jacobian  is  now  given  by 


[J]  = 


1   .   V  .T  .1  ^   1    d  P,  d^  +  /  G  ->-   ^    d  P3  d^  -   -   P,  d- 


2  A 


c,  b  P^  d  + 


4  A 


12 


4  A 


2  A 


-4  - 


^   dpj+   ^   bP^/  t.  ^l_b^  +  _i_dPg/  tf  c^b^ 


2  A 


2  A 


AA, 


4  A 


1  T    1 

d  Pc  c,  b^  +  _ 

"2"a"    ^   ^       2 


4  A 


2  b  t J  d  Pg  d"^  c^ 


2  V.T 


b  t^  d  P,  d^  t.  cf  b'  +  b  t'  d  P.  d^  t.  "1   b 


A  A, 


(4.68) 
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It  should  be  noted  that  the  Jacobian  is  symmetric.  The 
symmetry  is  obvious  in  terms  two,  three,  four,  seven,  nine,  eleven, 
and  twelve  of  equation  (4.68).  Terms  five  and  six,  and  terms  one  and 
nine  are  the  sum  of  a  matrix  and  its  transpose,  which  is  symmetric. 
Since  the  Jacobian  was  obtained  from  the  element  matrix  it  has  to  be 
calculated  along  with  the  element  properties  and  assembled. 

In  order  to  solve  the  problem  iteratively  the  residual  equation 

is  given  by 

[J]  5%  -  [RES],  (^-69) 

where 

4-^1-4  -  5tf.  (^-70) 

The  superscripts  in  equation  (4.70)   denote  the  iteration  numbers . 
Equation  (4.69)  can  now  be  rewritten  as 

[J]  4""^  "  ^^^  4  ■  ^^^^^-  ^""'^^^ 

The  assembled  form  of  equation  (4.71)  can  now  be  solved  iteratively. 
The  system  of  equations  in  each  iteration  were  solved  by  a  FORTRAN 
routine  SOLVE.  The  results  and  conclusions  are  presented  and  analysed 
in  Chapter  5.   Recommendations  for  further  study  are  also  presented. 
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CHAPTER  5 

RESULTS  AND  RECOMMENDATIONS 

In  this  work  a  continuum  approach  to  the  minimum  time  problem  has 
been  investigated  to  approximately  determine  the  isochrone 
distribution  in  the  phase  plane  of  a  double  integrator  problem.  The 
continuum  relations  were  derived  and  an  approximate  solution  was 
obtained  using  the  finite  element  method.  The  calculations  were  done 
by  FORTRAN  77  programs  implemented  on  an  IBM  370  VM/CMS  computer  and 
an  IBM  PC  microcomputer. 

A  four  by  four  square  region  was  selected  to  test  this  technique. 
A  number  of  different  grids  were  used  in  the  calculation  ranging  from 
very  coarse  to  very  fine.  Two  such  cases  are  presented  here.  Figure 
(5.1)  shows  the  isochrone  distribution  for  a  21  by  21  grid  while 
figure  (5.2)  shows  the  same  distribution  for  a  41  by  41  grid. 

The  abrupt  change  in  the  direction  of  the  isochrones  in  these  two 

figures  allows  one  to  approximately  trace  the  switching  curve.   The 

development  of  the  switching  line,  X2  ,  causes  the  switching  curve  to 

s 

have  a  tendency  to  prefer  a  horizontal  distribution.   The  true 

switching  curve  passes  through  the  lower  righthand  and  upper  lefthand 

corners  of  the  grid.   The  41  by  41  grid  shows  that  as  the  element 
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sizes   is   reduced  the  switching  curve  approaches  the  expected 

orientation. 

The  solutions  are  also  presented  in  tabular  form  for  four  cases 
along  with  the  actual  solutions.  A  comparison  of  the  actual  solution 
to  the  approximate  solution  shows  that  the  solutions  with  fine  grids 
give  very  good  results. 

Better  isochrone  distributions  are  possible  if  the  grid  is  made 
to  conform  to  the  switching  curve,  however,  this  defeats  the  purpose 
of  the  investigation. 

The  boundary  conditions  on  the  outer  grid  surface  were  varied  to 
study  the  performance.  The  only  set  of  boundary  conditions  which 
provided  a  reliable  solution  was  to  specify  the  true  final  time  at 
each  boundary  node.  That  this  is  a  necessary  condition  for  the 
solution  was  verified  by  deriving  the  differential  equations 
describing  the  variations  of  the  final  time  with  respect  to  the  state 
variables,   x^  and  X2 ,  from  equation  (3.31).   These  two  differential 

equations  are  separable  and  each  requires  two  boundary  conditions 
before  a  unique  solution  is  possible. 

An  interesting  condition  occurs  when  all  the  boundary  conditions 
are  removed  from  the  problem  except  the  boundary  condition  at  the 
origin.   The  solution  which  results  is 

t^ix^.x^)   -   -  UX2  (51) 

where  u  is  negative  for  positive  X2  and  positive  otherwise.   Equation 
(5.1)  produces  the  control  necessary  to  reduce  the  velocity  to  zero  in 
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minimum  time  regardless  of  the  initial  and  final  value  of  x^ .   An 
examination  of  equation  (5.1)   and  equation  (4.11)  shows  that  this 
simple   functional   expression  causes  equation   (4.11)   to  vanish 
entirely,  thus  rendering  equation  (4.11)  an  absolute  minimum. 
Conclusions 

Conclusions  reached  in  this  Investigation  are: 

1.  The  present  formulation  is  necessary  but  not  sufficient  to 
uniquely  specify  the  minimum  time  control;  other  information 
is  necessary  in  order  to  obtain  a  solution. 

2.  Coarser  grids  require  very  little  time  to  produce  an 
approximate  solution,  a  point  which  has  provided  motivation  to 
continue  the  work  in  this  area. 

3.  Linear  and  nonlinear  problems  present  the  same  level  of 
complexity  when  approached  from  this  point  of  view. 

4.  The  method  can  be  used  to  provide  closed  loop  minimum  time 
controls  by  determining  the  position  of  the  system  at  each 
instant  of  time  and  calculating  the  control  at  that  point. 

Problem  areas  requiring  additional  investigation  include: 

1.  The  determination  of  additional  continuum  relations  which 
would  provide  a  unique  solution  without  the  need  to  specify 
the  solution  on  the  boundary. 

2.  Developing  a  means  to  treat  systems  of  order  greater  than  two. 

3.  Determining  an  alternative  analysis  by  which  the  continuum 
condition  is  determined  directly  from  the  minimum  time 
functional. 
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4.  Testing  the  control  produced  by  this  or  similar  methods  to  see 
if  there  is  not  excessive  chattering. 
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(-2.-2) 


(2,2) 


Fig-re  5.1:  Isochrone  Discribucion  from  the  21x21  Finite 
Element  Grid 


(2,2) 


c-2.-:-. 


Figure  5.2:  Isochr-ne  Distribution  from  the  41x41  Finite 
Element  Grid 
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2.0  4.0  4.8  5.5  6.0 

2.2  1.4  2.4  3.4  4.2 

2.8  2,0  0.0  2.0  2.8 

4.2  3.4  2.4  1.4  2.2 

6.0  5.5  4.8  4.0  2.0 


Table  5.1a:  Approximate  Solution  from  the  Continuum  Approach 
For  the  5x5  Grid 


2.0  4.0  4.8  5.5  6.0 

2.2  1.9  2.2  3.4  4.2 

2,8  1.9  0.0  1.9  2.8 

4.2  3.4  2.2  1.9  2.2 

6.0  5.5  4.8  4.0  2.0 


Table  5.1b:  Actxial  Solution  from  the  Continuum  Approach 
For  the  5x5  Grid 
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2.0  3.3  3.9  4.3  4.7  5.0  5.3  5.5  5.8  6.0 

2.0  2.2  2.5  2.9  3.4  d.'9  4.3"'  4.6  4.9  5.1 

2.1  1.9  1.8  1.9  2.3  2.9  3.4  3.7  4.1  4.3 
2.3  2.0  1.6  1.3  1.3  1.9  2.5  3.0  3.3  3.6 
2.6  2.3  1.9  1.4  O.a  1.2  1.9  2.3  2.7  3.1 
3.1  2.7  2.3  1.9  1.2  0.8  1.4  1.9  2.3  2.6 
3.6  3.3  3.0  2.5  1.9  1.3  1.3  1.6  2.0  2.3 
4.3  4.1  3.7  3.4  2.9  2.3  1.9  1.8  1.9  2.1 
5.1  4.9  4.6  4.3  3.9  3.4  2.9  2.5  2.2  2.0 
6.0  5.8  5.5  5.3  5.0  4.7  4.3  3.9  3.3   2.0 


Table  5.2a:  Approximate  Solution  from  the  Continuum  Approach 
For  the  10x10  Grid 


2.0  3.3  3.9  4.3  4.7  5.0  5.3  5.5  5.8  6.0 

2.0  1.8  2.2  3.0  3.5  3.9  4.3  4.6  4.9  5.1 

2.1  1.8  1.5  1.2  2.4  2.9  3.4  3.7  4.1  4.3 
2.3  2.0  1.6  1.2  0.7  2.0  2.6  3.0  3.3  3.6 
2.6  2.3  1.9  1.4  0.8  1.2  1.9  2.4  2.7  3.1 
3.1  2.7  2.4  1.9  1.2  0.8  1.4  1.9  2.3  2.6 
3.6  3.3  3.0  2.6  2.0  0.7  1.2  1.6  2.0  2.3 
4.3  4.1  3.7  3.4  2.9  2.4  1.2  1.5  1.8  2.1 
5.1  4.9  4.6  4.3  3.9  3.5  3.0  2.2  1.8  2.0 
6.0  5.8  5.5  5.3  5.0   4.7  4.3  3.9  3.3  2.0 


Table   5.2b:    Actual   Solution   from   the   Continuum  .Approach 
For   the   10x10   Grid 
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APPENDIX  1 


CONSTRAINED  OPTTMIZATTOW  OF  A  FUNCTIONAL 

Consider  the  functional,  J,  given  by  equation  (2.18) 

J(x(t),  x(t).  u(t),  A(t).  t)-;^^i(x(t).  x(t).  u(t).  A(t).  t)dt  (Al.l) 

where  the  initial  time.  tg.  and  initial  states.  xCtp),  are  fixed  and 

the  final  time,   t^.   is  free  to  vary.   The  functional ,  J ,  can  be 


minimized  by  setting  the  variation  SJ   to  zero.  The  variation,  5 J ,  is 
given  by 
t, 


5J  =  Jt 


0 


^  5x  +!!!  5x  +  !!!  5u  +  !L  5A 
ax       du  d\ 

ax     ~ 


t^+dtp 
dt  +  ;^     g  dt.  (A1.2) 
^f 


Integrating  the  first  term  of  equation  (Al.2)  by  parts,  we  get 


-T 
5J  =  ff  (x(t).  x(t),  u(t),  A(t).  t) 

ax 


Sx(tf) 


+  g  (x(t),  x(t),  u(t),  A(t).  t) 


t-tf  ^^f 


+  //{[!!  (2i(t).  x(t).  u(t),  A(t).  t) 

^0  ax 
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t.   (—   (S(t),  x(t),  u(t),  A(t).  t))]  5x(t) 

dt   -■ 
3x 


+  !|  (x(t),  x(t),  u(t),  A(t).  t)  5u(t) 
dn 


+  !i  (x(t),  x(t),  u(t),  A(t).  t)  5A(t))  dt.     (A1.3) 
d\ 

Since  the  variation  5x  is  zero  at  the  initial  time  equation  (Al.3)  can 
be  rewritten  as 

,-T 


5J 


-  ^S  (x(t).  x(t),  u(t),  A(t),  t) 


d£ 


+  g  (x(t),  x(t),  u(t),  A(t).  t) 


5x(tj) 


t-tj  ^^f 


+  j/([!i  (s(t).  s(t),  u(t),  A(t),  t) 
'^o  ax 


^  (^  (x(t).  x(t),  u(t),  A(t).  t))]  5x(t) 

^t  ax 


+  fi  (x(t).  x(t),  u(t).  A(t),  t)  «u(t) 

au 


+  ?|  (x(t).  x(t),  u(t),  A(t),  t)  5A(t))  dt. 
dX 


(A1.4) 
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If  we  represent  the  quantity  dx(tj)   as   the  difference  between 
x(tj+dt^)  and  x(tj) .   Then  we  have 


or 


dx(tj)  -  5x(tj)  +  *(tj)  dt^ 
«S(tf)  -  dx(tj)  -  x(tj)  dtj 


(A1.5) 


Equation  (Al.2)  can  now  be  rewritten  as 

-T 
5J  -  ^   (x(t),   x(t),    u(t),   A(t),    t) 

ak 


t-t. 


dx(tp 


[!l  (x(t).   x(t),   u(t),   A(t),    t)   x(t)] 

ax 


dt. 


t-t. 


+  g   (x(t),    x(t),    u(t),   A(t).    t) 


t-t^  ^*^f 


+  //{[!!   (2i(t),    x(t),    u(t),    A(t),    t) 

0    ax 


i_   (fi   (x(t),    x(t),    u(t),   A(t),    t))]    Sx(t) 

d^    ax 


+  !f   (x(t),   x(t),    u(t).   A(t),    t)    5u(t) 

au 


+  !f   (x(t).    x(t),    u(t),   A(t),    t)    5A(t))    dt. 

dx 


(A1.6) 
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APPENDIX  2 

*    PROGRAM  MAIN  FOR  FINITE  ELEMENT  PACKAGE 

PARAMETER  (MBAND-43) 

PARAMETER  (IMAX-441) 

PARAMETER  (JMAX-IMAX*MBAND) 

PARAMETER  (NDMAX-6) 

PARAMETER  (NELMS-450) 

PARAMETER  (NTLEN-NELMS*NDMAX) 

PARAMETER  (NJ=NTLEN/3) 

PARAMETER  (IELMAX-900) 

REAL*4  X(IMAX)  ,Y(IMAX)  ,A(JMAX)  ,RSV(IMAX)  ,Em(NDMAX,NDMAX)  , 

+  PRSV(NDMAX).XI(3),YI(3),NT(NJ,3),ERROR,SUMTS, 

+  T ( IMAX) , TEMP , TOLD ( IMAX) , TTEST 

INTEGER*4  IDEBUG , NDS , NX , NY , NMAX/IMAX/ , MAXBND/MBAND/ , 

+  MAXMAT/JMAX/ , MAXEND/NDMAX/ , MAXELM/NELMS/ , 

+  MAXTL/NTLEN/ , NTABLE (NTLEN) , NUMEL , NELM , ITYPE , 

+  TWIDTH,NBANDW,ITABLE(NDMAX) ,NEQ, IT, IB(IELMAX) 

CHARACTER*!  ANS 

CHARACTER*80  FORM 

CHARACTER*2  NB 

EQUIVALENCE  (T(1),RSV(1)) 

COMMON  XL.YL.XH.YH 

READ(5,*)XL,YL,XH,YH 
10    IDEBUG-0 

WRITE(6,100)  ,      «„  „  xo  ,x 

100   FORMAT (LX,'  DO  YOU  WANT  THE  DEBUG  SWITCH  ON  (Y  OR  N  )?:') 

READ(5,110)ANS 
110   FORMAT (Al) 

IF(ANS.EQ.'Y')IDEBUG-1 

120   WRITE(6,130)  ^  ,^ 

130   FORMATdX,'  SPECIFY  THE  GRID  DENSITY  (NX  BY  NY  ) :  ) 
READ  (5,*)NX,NY 
IF(NX.LT.1.0R.NY.LT.1)THEN 

WRITE(6,140)  „  ,  ,^ 

140      FORMATdX,'  ERROR  -  GRID  PARAMETER  MUST  BE  1  OR  GREATER  /) 
GOTO  120 
ELSE  IF((NX+1)*(NY+1).GT.NMAX)THEN 

WRITE(6,150)  ^^„  ^^,^ 

150      FORMAT (IX,'  ERROR  -  REQUESTED  GRID  EXCEEDS  AVAILABLE  STORAGE') 

GOTO  120 
ENDIF 
CALL  NODE (X,Y,NMAX, NDS, NX, NY, IDEBUG) 

ITYPE-1 

CALL  GRIDLT (NX , NY , NTABLE , NELM , NUMEL , IDEBUG ) 
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END  IF 

CALL  BAND (NTABLE , TWIDTH , NUMEL , NBANDW , NELM , IDEBUG) 

NEQ-MAXMAT/NBANDW 

IF( (NBANDW . GT . MAXBND) . OR . (NEQ . LT . NDS )) THEN 

WRITE (6,*) 'ERROR  -  INSUFFICIENT  MEMORY' 
GOTO  10 
ENDIF 
TTEST-1 . 
DO  500  I-1,NDS 
T(I)-0. 
500   CONTINUE 

DO  900  IT-1,20 

IF(TTEST.GT.1.E-6)THEN 
DO  550  I-1,NDS 

T0LD(I)-T(I) 

^^°  CALL  BUILD (NTABLE, NUMEL, TWIDTH, NELM, A, RSV, NEQ. NBANDW, 

+  NDS , ITYPE , X , Y , NMAX , ELM , PRSV , MAXEND , IDEBUG , ITABLE , 

+  IB, TOLD, IT) 

CALL  BC (A, NBANDW, NEQ, RSV, NMAX, NDS, X.Y, IDEBUG) 

CALL  SOLVE (A , NBANDW , NEQ , RSV , NMAX , NDS ) 

ERROR-O . 0 

SUMTS-0 . 0 

DO  600  I-1,NDS 

ERROR-ERROR+ ( T ( I ) - TOLD ( I ) ) **2 

SUMTS-SUMTS+T ( I ) *T ( I ) 

600  CONTINUE 

TTEST-SQRT ( ERROR/SUMTS ) 

WRITE(6,610)IT,TTEST 
610  F0RMAT(1X,I5,4X,'TEST  -  ',E16.8) 

DO  700  J-1,NY+1 

WRITE(6,950)(T((J-1)*(NX+1)+I),I-1.NX+1) 

950  FORMAT(1X.24F5.1.10(/1X,24F5.1)) 

700  CONTINUE 

ENDIF 
900   CONTINUE 

NBA-NBANDW-1 
WRITE (NB, 25) NBA 
25    F0RMAT(I2) 

WRITE(6,*)'  THE  SOLUTION  IS  :' 
F0RM-'('//NB//'(lX,F4.1)/)' 
WRITE(6,F0RM)(RSV(I),I-1,NDS) 
DO  210  J-1,3 

DO  220  I-1,NJ 

NT(I,J)-NTABLE(NJ*(J-1)+I) 
220      CONTINUE 
210   CONTINUE 
ERROR-O . 
SUMTS-0 . 
DO  1100  I-1,NDS 

TEMP-GETTEE(X(I) ,Y(I) ) 
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ERROR-ERROR+ (TEMP - T ( I ) ) **2 
SUMTS-SUMTS+TEMP*TEMP 

1100  CONTINUE 

ERROR-SQRT ( ERROR/SUMTS ) 
WRITE (6, 1150) ERROR 
1150  F0RMAT(1X,'  ERROR-  ',E16.8) 

WRITE (6, 2000) 
2000  F0RMAT(1X,'D0  YOU  WANT  MAKE  A  PLOT  :') 
READ(5,110)ANS 
IF(ANS.EQ.'Y')THEN 
2050     WRITE(6,2100) 

2100     FORMAT (IX, 'HOW  MANY  CONTOURS  DO  YOU  WANT:') 
READ(5,*)I 
IF(I.LE.0)GOTO  2050 
NUMEL-MAXTL/3 

CALL  GRIDLT(NX.NY.NTABLE,NELM,NUMEL,IDEBUG)      ,,^,^^,^,0, 
CALL  PLT(X,Y,RSV,NMAX,NDS,NTABLE,NELM,NUMEL,A(l),A(MAXMAT/2) 

+  MAXMAT/ ( 2*1) , IDEBUG , I ) 

ENDIF 

WRITE (6, 10000)  ^^    ,^ 

10000  FORMAT (IX,'  DO  YOU  WANT  TO  RUN  ANOTHER  CASE   ) 
READ(5,110)ANS 
IF(ANS.EQ.'Y')THEN 

GOTO  10 
ENDIF 
STOP 
END 

REAL  FUNCTION  GETTEE(X,Y) 
REAL*4  X.Y.U 

U-1. 

IF((Y.GT.0.0  .AND.  X.GT. (-Y*Y*0. 5) ) .OR. 
+   (Y.LE.0.0   .AND.  X.GT. (Y*Y*0. 5) ) )U— 1. 
GETTEE— U*Y+2 . *SQRT ( -U*X+Y*Y*0 . 5 ) 
RETURN 
END 
*     SUBROUTINE  FOR  NUMBERING  NODES 

SUBROUTINE  NODE (X, Y, NMAX, NDS , NX, NY, IDEBUG) 
REAL*4  X(NMAX),Y(NMAX) 
INTEGER*4  NMAX, NDS , NX, NY, IDEBUG 
COMMON  XL,YL,XH,YH 
NDS-(NX+1)*(NY+1) 
DELTAX- (XH -XL) /FLOAT (NX) 
DELTAY- (YH-YL) /FLOAT (NY) 
IF (NY. GE. NX) THEN 
DO  10  I-1,NDS 

KK-M0D(I,(NX+1)) 
IF(KK.EQ.O)THEN 

X(I)-DELTAX*FLOAT(NX)+XL 

ELSE 
X(I)-DELTAX*FLOAT(KK- 1)+XL 

ENDIF 

78 


IY-(I-1)/(NX+1) 
Y(I)-YH-DELTAY*FLOAT(IY) 

10       CONTINUE 
ELSE 
DO  20  I-1,NDS 

IX-(I-1)/(NY+1) 

KK-M0D(I,(NY+1)) 

IF(KK.EQ.0)THEN 

Y(I)-YH-DELTAY*FLOAT(NY) 

ELSE 

Y ( I ) -YH - DELTAY*FL0AT ( KK - 1 ) 
ENDIF 

X ( I ) -DELTAX*FL0AT ( IX) +XL 
20       CONTINUE 
ENDIF 
IF(IDEBUG.NE.O)THEN 

WRITE(9,*)'THE  X  AND  Y  VALUES  ARE' 
DO  30  I-1,NDS 

WRITE(9,40)I,X(I),Y(I) 
40  FORMAT(4X,I4,2(4X,F10.6)/) 

30       CONTINUE 
ENDIF 
RETURN 
END 
*     SUBROUTINE  FOR  FORMING  THE  NODE  TABLE  FOR  EACH  ELEMENT 
SUBROUTINE  GRIDLT (NX , NY , NTABLE , NELM , NUMEL , IDEBUG) 
INTEGER*4  NX , NY , NTABLE (NUMEL , 3 ) , IDEBUG , Nl , N2 , N3 , N4 
NELM-0 
DO  10  I-1,NY 

DO  20  J-  1,NX 

IF (NX. GE. NY) THEN 

N1-(I-1)*(NX+1)+J 
N2-N1+1 
N4-I*(NX+1)+J 
N3-N4+1 
ELSE 

N1-(J-1)*(NY+1)+I 
N4-N1+1 
N2-J*(NY+1)+I 
N3-N2+1 
ENDIF 

NELM-NELM+1 

IF (NELM. GT. NUMEL) GOTO  50 
NTABLE ( NELM, 1)-N1 
NTABLE (NELM, 2 )-N4 
NTABLE (NELM, 3 )-N2 
NELM-NELM+1 

IF (NELM. GT. NUMEL) GOTO  50 
NTABLE (NELM, 1)-N2 
NTABLE (NELM, 2 )-N4 
NTABLE (NELM, 3 )-N3 
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20       CONTINUE 
10    CONTINUE 

TFfTDEBUG  NE  O^THEN 

WRITE(9.*)'  THE  NUMBER  OF  ELEMENTS  REQD.  FOR  THE  PROBLEM  :', 

+  '   NELM 

WRITE(9,25) 
25       FORMAT(//'  THE  NODE  TABLE  IS  :'//) 

WRITE(9  ,  30)1  .NTABLEd .  1)  .NTABLECI ,  2)  ,NTABLE(1 .  3) 
30  F0RMAT(1X,I4,3(2X,I4)/) 

40       CONTINUE 

ENDIF 

RFTURN 
50    WRITE(5,*) 'ERROR  ***  NUMBER  OF  ELEMENTS  REQD.  EXCEEDS  MEMORY 

+       '    RESERVED' 

STOP 

END 
*     SUBROUTINE  FOR  BUILDING  THE  GLOBAL  MATRIX 

SUBROUTINE  BUILD(NTABLE,NUMEL,TWIDTH,NELM.A,RSV.NEQ,NBANDW  NDS , 
+  ITYPE,X,Y.NMAX,ELM,PRSV,MAXEND,IDEBUG,ITABLE, 

TR  TOT  n  TT^ 
^INTEGER*4  TWIDTH, ITABLE(MAXEND) ,NTABLE(NUMEL,TWIDTH) . IT , IB(NUMEL) 

REAL*4  ELM(MAXEND,MAXEND) ,PRSV(MAXEND) ,RSV(NMAX) ,A(NEQ,NBANDW) , 
+       X(NMAX),Y(NMAX),TOLD(NMAX) 
DO  10  I-1,NDS 

DO  20  J-1,NBANDW 
A(I.J)-0. 
20       CONTINUE 
RSV(I)-0. 
10     CONTINUE 

DO  30  I-1,NELM 

DO  40  J-1, TWIDTH 

ITABLE ( J ) -NTABLE ( I , J ) 
40       CONTINUE 

Xl-X(  ITABLE  (D) 
X2-X(ITABLE(2)) 
X3-X(ITABLE(3)) 
Yl-Y ( ITABLE (1)) 
Y2-Y(ITABLE(2)) 

SlL^UNTRI(Xi!x2.X3.Y1.Y2,Y3,ELM,PRSV,MAXEND,IDEBUG  IB(Ih 
+  T0LD(ITABLE(1)).T0LD(ITABLE(2)),T0LD(ITABLE(3)),IT) 

DO  60  K-1, TWIDTH 
IJ-ITABLE(K) 
DO  50  J-1, TWIDTH 

II-ITABLE(J) 

IF(II.LT.IJ)GOTO  50 

IK-II-IJ+1 
A(IJ,IK)-A(IJ,IK)+ELM(K,J) 

50  CONTINUE 

RSV ( IJ ) -RSV ( I J ) +PRSV (K) 
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60 
30 


21 
22 


C 

C 
C 
C 
C 
C 
C 
C 


CONTINUE 
CONTINUE 
IF(IDEBUG.NE.O)THEN 

WRITE(9,*)'  THE  GLOBAL  MATRIX  IS  :' 

WRITE(9,21)((A(I,J),J-1.TWIDTH).I-1,NDS) 

F0RMAT(3(F8.4,1X)/) 

WRITE(9,22)(RSV(I),I-1,NDS) 

F0RMAT(5(F8.4)/) 
ENDIF 
RETURN 
END 

SUBROUTINE  LINTRI (XI , X2 , X3 , Yl , Y2 , Y3 , ELM , ERSV , NM , IDEBUG , IB , 
&  T1,T2,T3,IT) 

LINTRI  FOR  IS  FOR  THE  SOLUTION  OF  THE  TRANSVERSALITY  CONDITION 
IN  THE  XI -X2  PHASE  PLANE  FOR  THE  DOUBLE  INTEGRATOR  PROBLEM.   THE 
EQUATIONS  ARE  NONLINEAR  AND  THE  SOLUTION  IS  OBTAINED  THROUGH 
NEWTON  METHOD   THIS  SUBROUTINE  BUILDS  THE  ELEMENT  JACOBIAN  MATRIX 
AND  EVALUATES  THE  RESIDUALS  ON  AN  ELEMENT  BASIS.   INCLUDED  WITH  THE 
THE  TRANSVERSALITY  CONDITION  IS  A  CLOSED  FORM  EXPRESSION  FOR  THE 
MINIMUM  TIME  FUNCTIONAL.   THE  ELEMENT  EQUATIONS  ARE  DETERMINED  FROM 
THE  LEAST  SQUARES  PROCESS. 


INTEGER*4  NM,  IDEBUG,  IB,  IT 

REAL*4  XI,  X2,  X3 ,  Yl,  Y2,  Y3,  ELM(NM,NM),  ERSV(NM) 

REAL*4  Tl,  T2,  T3 

LOCAL  VARIABLES 

INTEGER*4  I,  J,  lU,  lUl 

REAL*4  A(3),  B(3) ,  C(3) ,  U(3) .  X(3),  Y(3),  LI,  L2 ,  DEL2 ,  X2S 

REAL*4  UJ,  Al,  UINT,  X2UINT,  DUINT,  DDUINT,  DX2U,  DEL,  DDX2U 

REAL*4  Hll,  H12,  H22,  CI,  C2,  C3 ,  C4,  C5 ,  C6 ,  C7 ,  C8 ,  09 

REAL*4  T(3),  XAVE,  YAVE,  H(3,3),  Gil,  G12,  G22,  F(3),  XXX 


A(l) 
A(2) 
A(3) 
B(l) 
B(2) 
B(3) 
0(1) 
C(2) 
0(3) 
Y(l) 
Y(2) 
Y(3) 
X(l) 
X(2) 
X(3) 
T(l) 
T(2) 


X2*Y3  -  X3*Y2 
X3*Y1  -  X1*Y3 
X1*Y2  -  X2*Y1 
Y2  -  Y3 
Y3  -  Yl 
Yl  -  Y2 
X3  -  X2 
XI  -  X3 
XI 


-  X2 

-  Yl 

-  Y2 

-  Y3 

-  XI 

-  X2 

-  X3 

-  Tl 

-  T2 
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T(3)  -  T3 

XAVE  -  (Xl+X2+X3)/3.0 

YAVE  -  (Yl+Y2+Y3)/3.0 

DEL2  -  A(l)  +  A(2)  +  A(3) 

IF(DEL2  .LE.  0.0)  THEN 

WRITE(*,80)  DEL2  ,  ^^ ^  ^^ 

FORMAT (IX, 'ERROR  IN  GRID  :  2  *  TRIANGLE  AREA  -',F12.6) 

STOP 
ENDIF 


DEL  -  0.5  *  DEL2 

LI  -  (-B(1)*T1-B(2)*T2-B(3)*T3)/DEL2 

L2  -  (-C(1)*T1-C(2)*T2-C(3)*T3)/DEL2 

Hll  -  DEL2* ( Y1*Y1+Y2*Y2+Y3*Y3+Y1*Y2+Y1*Y3+Y2*Y3 ) *2 . 0/24 . 0 

H22  -  DEL 
DO  100  I  -  1,  3 
U(I)  -  1.0 

IF(  IT  .EQ.  1)  THEN  „  ^^   ^„ 

IF(  YAVE  GT.  0.0  .AND.  XAVE  .GT.  (-YAVE*YAVE*0. 5)  .OR. 
&  YAVE  .LE.  0.0  .AND.  XAVE  .GT.  (  YAVE*YAVE*0 . 5 ) ) THEN 

U(I)  -  -1.0 
ENDIF 
ELSE 

IF(  L2*(1.0  -  L1*Y(I))  .LT.  0.0  )  U(I)  -  -1.0 

ENDIF 
100   CONTINUE 

DO  200  I  -  1,  3 

DO  180  J  -  1,  3 

H(I,J)  -  2.0  *  X(I)  *  B(J)  +  Y(I)  *  C(J) 
180       CONTINUE 
200   CONTINUE 

DO  220  I  -  1,  3 
F(I)  -  0.0 
DO  210  J  -  1,  3 

F(I)  -  F(I)  +  H(J,I) 
210       CONTINUE 

220   CONTINUE 

Gil  -  DEL2*(Xl*Xl+X2*X2+X3*X3+Xl*X2+Xl*X3+X2*X3)/3.0 
G12  -  DEL2*(2.0*(X1*Y1+X2*Y2+X3*Y3)+X1*Y2+X2*Y1+X1*Y3+X3*Y1+ 
^  X2*Y3+X3*Y2)/12.0 

G22  -  DEL2* ( Y1*Y1+Y2*Y2+Y3*Y3+Y1*Y2+Y1*Y3+Y2*Y3 ) /12 . 0 

DO  250  I  -  1,  3 

ELM(iTj)'-  DEL2/24.0  -  (H(I,J)+H(J,I)+F(I)+F(J))/24  0 
ELM(I  J)  -  ELM(I.J)  +  (B(I)*B(J)*G11+(B(I)*C(J)+B(J)*C(I)) 
4  *G12+C(I)*C(J)*G22)/(DEL2*DEL2) 

IF(I  .EQ.  J)  THEN 

ELMd.I)  -  ELMd.I)  +  DEL2/24.0 
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ELSE 

ELM(J,I)  -  ELM(I,J) 
ENDIF 
240       CONTINUE 
250   CONTINUE 

DO  270  I  -1,  3 

ERSV(I)  -  0.0 
DO  260  J  -  1,  3 

ERSV(I)  -  ERSV(I)  -  ELM(I,J)*T(J) 

260       CONTINUE 
270   CONTINUE 

IF(  U(1)*U(2)  .LT.  0.0  .AND.  IT  .NE.  1  .OR. 
&  U(1)*U(3)  .LT.  0.0  .AND.  IT  .NE.  1  .OR. 
&    U(2)*U(3)  .LT.  0.0  .AND.  IT  .NE.  1  )   THEN 

IF(  U(l)  .NE.  U(2)  .AND.  U(l)  .NE.  U(3)  )  lU  -  1 

IF(  U(2)  .NE.  U(l)  .AND.  U(2)  .NE.  U(3)  )  lU  -  2 

IF(  U(3)  .NE.  U(l)  .AND.  U(3)  .NE.  U(2)  )  lU  =  3 

X2S  -  1.0  /  LI 

UJ  -  U(IU) 

lUl  -  lU  -  1 

IF(IU1  .EQ.  0)  lUl  -  3 

Al     -  DEL  *  (X2S  -  Y(IU))/(Y(IU1)  -  Y(IU)) 

UINT   -  UJ  *  (2.0*A1  -  DEL) 

X2UINT  -  UJ*(2.0*Al*(2.0*X2S+Y(IU))-DEL*(Yl+Y2+Y3))/3.0 

DUINT  -  UJ/((Y(IU1)-Y(IU))*L1*L1) 

DDUINT  -  DUINT  /  (DEL  *  LI)  ..  „  «^.-,x  / 

DX2U   -  UJ*(DEL*(2.0*X2S+Y(IU))/(Y(IU1)-Y(IU))+2.0*A1)/ 

&  (3.0*DEL*L1*L1) 

DDX2U  -  2.0*UJ/(3.0*(Y(IU1)-Y(IU))*DEL*(L1**4))+DX2U/(DEL*L1) 

CI  -  (DUINT+X2UINT/DEL2)/DEL2 

C2  -  1.0/(DEL2*DEL2) 

C3  -  -DX2U/DEL2 

C4  -  -L2*DDUINT  +  L1*L2*DDX2U 

C5  -  CI  +  LI  *  C3 

C6  -  C2*H11  +  C3*2.0*L2  +  C4 

P7  ^  P9^H22 

C8  -  ( Y1+Y2+Y3 ) /6 . 0 - L1*H11/DEL2+ (Ll*DX2U-X2UINT/DEL2 - DUINT ) *L2 

C9  -  (UINT  -  L2*H22  -  Ll*X2UINT)/DEL2 

DO  500  I  -  1,  3 

ERSV(I)  -  ERSV(I)  -  C8*B(I)-C9*C(I) 

DO  300  J  -  I,  3 

ELM(I,J)-ELM(I.J)+C5*(C(I)*B(J)+C(J)*B(I))+ 
&  C6*B(I)*B(J)+C7*C(I)*C(J) 

ELM(J,I)-ELM(I,J) 
300  CONTINUE 

DO  400  J  -  1,  3 

ERSV(I)  -  ERSV(I)  +  ELM(I,J)*T(J) 

400  CONTINUE 
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500       CONTINUE 
ELSE 

IF(  IT  .EQ.  1)  THEN 

IF("yAVE  GT  0.0  .AND.  XAVE  .GT.  (-YAVE*YAVE*0 . 5)  .OR. 
&  YAVE  '.LE.  0.0  .AND.  XAVE  .GT.  (  YAVE*YAVE*0  .  5 ) ) THEN 

UJ  -  -1.0 
ENDIF 
ELSE 

UJ  -  U(l) 
ENDIF 
H12  -  UJ*DEL*(Yl+Y2+Y3)/3.0 

UINT  -  UJ*DEL 

^°  ^ERSV(I)^-  ERSV(I)  -  B(I)*(Yl+Y2+Y3)/6.0  -  G(I)*UINT/DEL2 

DO  600  J  -  1,  3 

XXX  -  (B(I)*B(J)*H11+(C(I)*B(J)+C(J)*B(I))*H12+ 

&  C(I)*C(J)*H22)/(DEL2*DEL2) 

ERSV(I)  -  ERSV(I)  -  XXX  *  T(J) 
ELM(I,J)  -  ELM(I,J)  +  XXX 
600  CONTINUE 

700       CONTINUE 

DO  900  I  -  1,  3 

DO  800  J  -  1,  3 

ERSV(I)  -  ERSV(I)  +  ELM(I,J)  *  T(J) 
800  CONTINUE 

900       CONTINUE 
ENDIF 

IF(IDEBUG  .  NE.  0)  THEN 

WRITE(*,2000)  XI,  X2,  X3 ,  Yl,  Y2,  Y3 
2000      FORMAT (IX, 'LINEAR  TRIANGLE' , 6F12 . 5) 
DO  2500  I  -  1,  3 

WRITE(*,2200)  (ELM(I,J),J-1,3),  ERSV(I) 
2200         F0RMAT(1X,3F12.5,  lOX,  F12.5) 
2500      CONTINUE 
ENDIF 

RETURN 
END 
*     SUBROUTINE  FOR  INCLUDING  THE  BOUNDARY  CONDITIONS 

SUBROUTINE  BC ( A , NBANDW , NEQ , RSV , NMAX , NDS , X , Y , IDEBUG ) 

INTEGERS  NBANDW , NEQ , NMAX , NDS , IDEBUG 

REAL*4  A(NEQ, NBANDW) ,RSV(NMAX) ,X(NDS) ,Y(NDS) 

CHARACTER*80  FORM 

CHARACTER*2  NB 

COMMON  XL,YL,XH,YH 

DO  10  I-1,NDS  ,  ^  ^^  „„ 

IF(((ABS(X(I)-XL)).LE.l.E-3).OR.((ABS(X(I)-XH)).LE.l.E-3).OR. 

+    ((ABS(Y(I)-YL)).LE.l.E-3).OR.((ABS(Y(I)-YH)).LE.l.E-3).OR. 

+    ((ABS(X(I)).LE.1.E-3).AND.((ABS(Y(I))).LE.1.E-3)))THEN 
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IF(((Y(I).GT.0.0).AND.(X(I).GT.(-Y(I)*Y(I)*0.5))).OR. 
+  ((Y(I).LE.  0.0).AND.(X(I).GT.(Y(I)*Y(I)*0.5))))THEN 

U— 1.0 
ENDIF 

SQA— U*X(I)+Y(I)*Y(I)*0.5 
IF(SQA.LE.l.E-6)SQA-0. 
RSV(I)— U*Y(I)+2.*SQRT(SQA) 
IF(RSV(I).LT.O.)RSV(I)— RSV(I) 
RSV(I)-RSV(I)*A(I,1)*1.E20 
A(I,1)-A(I,1)*1.E20 
ENDIF 
10    CONTINUE 

WRITE(NB,5)NBANDW 
5     FORMAT (12) 

FORM-' ( ' //NB//' (E12 . 3 , IX)/) ' 
IF(IDEBUG.NE.O)THEN 

WRITE(9,*)'  THE  GLOBAL  MATRIX  IS  :' 
WRITE(9 , FORM) ( (A(I , J) , J-1 .NBANDW) , I-1,NDS) 

WRITE(9,112)(RSV(I).I-1,NDS) 
112      F0RMAT(5(E14.4,1X)/) 

ENDIF 

RETURN 

END 
*     SUBROUTINE  FOR  SOLVING  THE  GLOBAL  MATRIX 

SUBROUTINE  SOLVE  (A,  HBANDW,  NEQ,  RSV,  NMAX,  NDS) 

INTEGER*4  HBANDW,  NEQ,  NMAX,  NDS,  I,  J,  K,  MAXCOL 
REAL*4     A(NEQ, HBANDW),  RSV (NMAX) ,  DENOM,  TERM,  SUM 
C  REDUCE  MATRIX 
C   LOOP  ON  COLUMNS  TO  BE  REDUCED 
DO  300  I  -  1,  NDS-1 

DENOM  -  1.0/A(I,1) 

MAXCOL  -  NDS  -  I  +  1 

IF (MAXCOL  .GT.  HBANDW)  THEN 

MAXCOL  -  HBANDW 
ENDIF 
DO  200  J  -  2,  MAXCOL 

IF(ABS(A(I,J))  .GT.  1.0E-20)  THEN 
TERM  -  A(I,J)  *  DENOM 
IF(ABS(TERM)  .GT.  l.OE-25)  THEN 
DO  100  K  -  J,  MAXCOL 

IF(ABS(A(I,K))  .GT.  l.OE-25)  THEN 

A(I+J-1,K-J+1)-A(I+J-1,K-J+1)-TERM*A(I,K) 

ENDIF 
100  CONTINUE 

RSV(I+J-1)  -  RSV(I+J-1)  -  TERM  *  RSV(I) 
ENDIF 
ENDIF 
200       CONTINUE 
300   CONTINUE 
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C   BACK  SUBSTITUTE 

DO  500  I  -  NDS,  1,  -1 
SUM  -  RSV(I) 
DO  400  J  -  2,  HBANDW 
K  -  I  +  J  -  1 
IF(K  .LE.  NDS)  THEN 

IF(ABS(RSV(K))  .GT.  l.OE-25)  THEN 

IF(ABS(A(I,J))  .GT.  l.OE-25)  THEN 

SUM  -  SUM  -  RSV(K)  *  A(I,J) 
ENDIF 
ENDIF 
ENDIF 
400       CONTINUE 

RSV(I)  -  SUM  /  A(I,1) 
500   CONTINUE 
RETURN 

*     SUBROUTINE  FOR  DIVIDING  THE  REGION  INTO  CONTOURS 

SUBROUTINE  PLT(X,  Y,  T,  NMAX,  NDS,  NTABLE,  NELM,  NUMEL,  XS , 
&  YS,  MAXP,  IDEBUG,  NLINE) 

INTEGERS      I.  NLINE.  NMAX,  NPP(IOO) ,  J.  K.  Kl,  K2 ,  IMIN   IMAX 
INTEGER*4      NDS,  IDEBUG,  NELM,  NUMEL,  MAXP.  NTABLE (NUMEL  3) 
REAL*4         X(NMAX),  Y(NMAX) ,  XS (MAXP, NLINE ) ,  YS (MAXP, NLINE) 
1^1*1  T(NMAX);  TMAX.  TOIN.  XX.  YY.  RATIO,  TINC.  ETMAX 

REAL*4         ETMIN 
CHARACTER*1    ANS 
IF (IDEBUG  .NE.  0)  THEN 
DO  10  I  -  1,  NELM 

WRITE(9,5)I.(NTABLE(I.J),J-1,3) 

5  F0RMAT(1X.I4.5X,3I6) 

10  CONTINUE 

DO  20  I  -  1,  NDS 

WRITE(9,15)I,X(I).Y(I).  T(I) 
15  F0RMAT(1X,I5,3F12.5) 

20        CONTINUE 
ENDIF 

11  TMAX  -  T(l) 
TMIN  -  T(l) 

DO  50  I  -  2,  NDS 

IF(T(I)  .GT.  TMAX)  TMAX  -  T(I) 

IF(T(I)  .LT.  TMIN)  TMIN  -  T(I) 

50    CONTINUE 

DO  1000  I  -  1,  NLINE 

NPP(I)  -  0 

TINC  -  FLOAT(2*I-l)*(TMAX-TMIN)/FLOAT(2*NLINE)  +  TMIN 

IF(IDEBUG  .NE.  0)  WRITE(9,55)  TINC 
55        FORMAT (IX, 'TINC  -  ',E16.8) 
DO  500  J  -  1,  NELM 

ETMAX  -  T(NTABLE(J,1)) 

ETMIN  -  T (NTABLE ( J, D) 
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IMAX  -  1 
IMIN  -  1 
DO  100  K  -  2,  3 

IF(T(NTABLE(J,K)).GT.  ETMAX)  THEN 
ETMAX  -  T(NTABLE(J,K)) 
IMAX  -  K 
ENDIF 

IF(T(NTABLE(J,K))  .LT.  ETMIN)  THEN 
ETMIN  -  T(NTABLE(J,K)) 
IMIN  -  K 
ENDIF 
100  CONTINUE 

IF(TINC  .GE.  ETMIN  .AND.  TING  .LE.  ETMAX)  THEN 

DO  400  K  -  1,  3 
Kl  -  K 
K2  -  K  +  1 
IF(K2  .EQ.  4)  K2  -  1 

IF(ABS(T(NTABLE(J ,K1) ) -T(NTABLE(J ,K2) ) ) 
&  LT  1  OE-5  *  ABS(TINC)  .AND.  ABS(T(NTABLE(J ,K1) ) - 

&  TING)  .LT.  l.OE-5  *  ABS(TING))  THEN 

XS(NPP(I)+1,I)  -  X(NTABLE(J,K1)) 
XS(NPP(I)+2,I)  -  X(NTABLE(J.K2)) 
YS(NPP(I)+1.I)  -  Y(NTABLE(J.K1)) 
YS(NPP(I)+2,I)  -  Y(NTABLE(J.K2)) 
NPP(I)  -  NPP(I)  +  2 
ELSE 

ETMIN  -  T(NTABLE(J,K1)) 
ETMAX  -  T(NTABLE(J,K2)) 
IF (ETMIN  .GT.  ETMAX)  THEN 
XX  -  ETMIN 
ETMIN  -  ETMAX 
ETMAX  -  XX 
ENDIF 

IF(TING  .GT.  ETMIN  .AND.  TING  .LT.  ETMAX) THEN 
RATIO  -  (TING-T(NTABLE(J,Kl)))/( 
&  T(NTABLE(J,K2))  -  T(NTABLE(J ,K1) ) ) 

XX  -  RATIO  *  (X(NTABLE(J,K2))- 
&  X(NTABLE(J.K1)))  +  X(NTABLE(J ,K1) ) 

YY  -  RATIO  *  (Y(NTABLE(J,K2))- 
&  Y(NTABLE(J.K1)))  +  Y(NTABLE(J ,K1) ) 

NPP(I)  -  NPP(I)  +  1 
XS(NPP(I),I)  -  XX 
YS(NPP(I),I)  -  YY 
ENDIF 
ENDIF 
400  CONTINUE 

ENDIF 
500       CONTINUE 

IF(NPP(I)  .GT.  MAXP)  THEN 
WRITE(6,600) 
600  FORMAT (IX, 'MORE  DATA  POINTS  REQUIRED  FOR  PLOT  THAN' , 
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&        •  AVAILABLE  -  USE  FEWER  LINES') 
READ(5,111)ANS 
GOTO  11 
ENDIF 
1000  CONTINUE 

CALL  PLOT(XS,YS,NPP,NMAX,T,NLINE,MAXP) 

112   WRITE(5,*)'  DO  YOU  WANT  ANOTHER  PLOT:' 

READ(5,111)ANS 
111   FORMAT (Al) 

IF(ANS.EQ.'Y')THEN 

WRITE (5,*)'  ENTER  NO.  OF  CONTOURS:' 

READ(5,*)NLINE 

GOTO  11 

ELSE  IF(ANS.NE.'N')THEN 

GOTO  112 
ENDIF 
RETURN 
END 

SUBROUTINE  T04014 
INTEGER  ARRAY(2)/27,49/ 
C      CHARACTER  ICHRS 
C     ICHRS  -  '  ' 

CALL  KAS 2 AM (2, ARRAY, ICHRS) 
WRITE(6,555)  ICHRS 
555  FORMATC  '  ,  A2) 
RETURN 
END 

SUBROUTINE  TOANSI 
INTEGER  ARRAY(2)/27,50/ 
C      CHARACTER  ICHRS 
C     ICHRS  -  '  ' 

CALL  KAS 2 AM (2, ARRAY, ICHRS) 
WRITE(6,555)  ICHRS 
555  FORMATC  '  ,  A2) 
RETURN 
END 
*     SUBROUTINE  FOR  PLOTTING  THE  ISOCHRONES 

SUBROUTINE  PLOT  (XS ,YS ,NPP,NMAX,T,NLINE,MAXP) 

REAL*4  XS(MAXP,NLINE) ,YS(MAXP,NLINE) ,T(NMAX) 

INTEGER*4  NLINE,NPP(100) 

COMMON  XL,YL,XH,YH 

CALL  T04014 

CALL  GRSTRT(4014,1) 

CALL  NEWPAG 

CALL  WINDOW ( XL, XH.YL.YH) 

CALL  VWPORT(0. ,130. ,0. ,100.) 

CALL  CMOPEN 

CALL  MOVE (XL, YL) 

CALL  CMCLOS 
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CALL  CMOPEN 

DO  10  I-1,NLINE 

DO  20  J-l,NPP(I)/2 

CALL  MOVE(XS(2*J-l,I),YS(2*J-l.I)) 
CALL  DRAW(XS(2*J,I).YS(2*J,I)) 
20       CONTINUE 
10    CONTINUE 

CALL  MOVE(XL,YL) 

CALL  DRAW(XH,YL) 

CALL  DRAW(XH,YH) 

CALL  DRAW(XL,YH) 

CALL  DRAW(XL,YL) 

CALL  CMCLOS 

CALL  GRSTOP 

CALL  TOANSI 

RETURN 

END 
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ABSTRACT 


A  continuum  approach  to  the  minimum  time  problem  has  been  taken 
in  order  to  approximately  determine  the  isochrone  distribution  in  the 
phase  plane  of  a  double  integrator  problem.  The  approximate  solution 
and  the  actual  solution  have  been  provided. 

In  this  analysis  the  system  state  variabiles  have  been  treated  as 
independent  variables  and  time  has  been  eliminated  from  the  analysis. 
This  choice  of  variables  enables  us  to  treat  both  linear  and  nonlinear 
systems  with  the  same  methods  of  solution.  Moreover,  the  costate 
variables  can  be  determined  directly  from  the  gradient  with  respect  to 
the  state  variables  of  the  final  time  at  any  point  in  state  space  in 
accordance  with  the  Hamilton-Jacobi-Bellman  equation.  The  continuum 
relations  have  been  derived  and  an  approximate  solution  has  been 
obtained  for  a  double  integrator  problem  by  using  a  finite  element 

technique . 

The  approximate  solution  has  been  compared  with  the  actual 
solution  and  a  plot  of  the  isochrones  has  been  provided. 


